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SUMMARY 


The numerical computation of unsteady airloads acting upon thin airfoils 
with multiple leading and trailing edge controls in two dimensional ventilated 
subsonic wind tunnels is studied. The foundation of the computational method 
is strengthened with a new and more powerful mathematical existence and 
convergence theory for solving Cauchy singular integral equations of the first 
kind, and the method of convergence acceleration by extrapolation to the limit 
is introduced to analyze airfoils with flaps. New results are presented for 
steady and unsteady flow, including the effect of acoustic resonance between 
ventilated wind tunnel walls and airfoils with oscillating flaps. 
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§li Introduction ^ 

The theory of aerodynamic interference is necessary to extrapolate wind tunnel 
test data to free flight conditions. Analytically the wind tunnel problem is more 
difficult them the free air problem because of the presence of the wind tunnel walls, 
although computationally it may be checked against known free air results. However, 
the physical validity of wind tunnel interference theories must be established by 
the more stringent comparison with experimental data for both the free air and wind 
tunnel conditions. 

Physically the nature of wind tunnel flow is highly complex so that consider- 
able simplification of the theory is necessary in order to produce useful results. 

Most existing wind tunnel theories, see e.g., [1], [2], [3], [4]^, are based on far 

field effects under the assumption of inviscid linear potential flow. Clearly vis- 
cosity plays a major role in boundary layer and wall effects, as do nonlinearities 
in transonic flow. Although significant progress has been made in the general field 
of computational fluid mechanics in the past decade, the inclusion of viscosity and 
nonlinearities entails an order of magnitude more computational expense than the 
linear inviscid theories. 

In shock free flow over moderately thin airfoils, good correlation with experi- 
ment ccin be obtained using linear potential flow theory, and there remain open signi- 
ficant research problems of practical importance. Among these are the proper form 
or forms of the boundary condition at ventilated tunnel walls, the efficient compu- 
tation of unsteady airloads especially when flaps are present, and effects of acoustic 
resonance between oscillating airfoils and ventilated wind tunnel walls. 

In this report we discuss the problem of predicting unsteady airloads on thin 
airfoils in two dimensional subsonic flow through ventilated wind tunnels. Our pre- 
vious work [5] , [6] is extended to include the problem of multiple leading and trail- 

ing edge controls, and to permit porous wind tunnel walls in steady flow. This work 
is based on the method of orthogonal polynomial pairs discovered originally by N.I. 
Akheizer [7] in the Soviet Union in 1945 and developed independently but later by 
S.R. Bland cind his coworkers for the practical solution of airfoil problems [8], [9], 

^The authors wish to acknowledge the valuable help of Dr. Sanford Davis, Ms. Theda 
Grinnell and Mesrs. Charles Doughty, Paul Kriner and Steven Sedlacek. 

^Numbers in square brackets refer to the bibliography found in the REFERENCES in 
order of first citation, and may also give the page, section or chapter of interest. 
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[10] . In section 2 below we set down the underlying partial differential equations 
and boundairy conditions- Section 3 describes Bland’s kernel of the integral equa- 
tion relating downwash to pressure for unsteady flow in slotted wall tunnels and 
section 4 describes the kernel for porous wall tunnels in steady flow. Sections 5 
and 6 review the theory of airfoil polynomials and Bland's collocation method. Al- 
though this is slightly repetitious of earlier work, we have simplified the notation 
somewhat and unified certain other concepts. Sections 7 and 8 discuss the calcula- 
tion of airloads and the method of representation of airfoil profiles possessing mul- 
tiple leading and trailing edge controls. Section 9 and the APPENDIX present instruc 
tions for the use of the computer program TWODI and a sample interactive input/output 
scenario. Sections 10 and 11 present a comprehensive analysis of convergence of col- 
location and other computational methods and of the method of extrapolation used to 
accelerate numerical convergence in the case of flaps. Section 12 shows how computa- 
tional efficiency can be further improved by identifying the singularities in the ker 
nel and integrating them separately. Sections 13, 14 and the APPENDIX present numer- 
ical results for steady flow in porous wall tunnels and for steady and unsteady flow 
over airfoils with flaps in ventilated tunnels for frequency ranges which include 
resonance between the oscillating airfoil and the wind tunnel walls. In section 15, 
we offer our conclusions based on the present study and recommendations for future 
research. 
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§2. Basic equations 

Consider a thin nearly planar airfoil undergoing simple harmonic motion at 
frequency qj rad/sec about the center plane of a two dimensional ventilated wind 
tiinnel (Figure 1). The flow is considered to be inviscid and strictly subsonic, 
and to be a small disturbance from the free stream flow. 

We point out that free air conditions are included as an important special 
case upon taking the walls to be infinitely far apart . 




Figure 1. Coordinate system and sign conventions 

In this study we are particularly interested in unsteady interference effects 
on airfoils with flaps but, in, general, we may allow the deflection amplitude func- 
tion h to represent upper, lower or mean airfoil profiles as well as arbitrary aero- 
elastic shapes, including rigid body and control surface deflections. 

Deflection is measured positive downward so that positive streamwise deriva- 
tives correspond to positive angles of attack. Lift will be measured positive up- 
ward, and pitching moment will be measured about the quarter chord, positive in the 
direction of increasing angle of attack (nose up) . 
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Except where specifically stated to the contrary, all quantities are dimension- 
less. Lengths are nondimensionalized by the airfoil semichord b, pressures by the 

1 P 

free stream dynamic pressure — p and, with the notable exception of Mach nvimber, 

2 00 OO 

velocities are nondimensionalized by the free stream velocity v^. Following cus- 
tomary practice, streamwise coordinates are measured positive downstream, having 
value -1 at the leading edge and +1 at the trailing edge. 

Denote by w the downwash amplitude nondimensionalized by the free stream veloc- 
ity, 

w(x) = + ik)h(x), (2.1) 

dx 

where k=o)b/v^ is the reduced frequency based on semichord. Denote by c() the pertur- 
bation velocity potential nondimensionalized by bv and denote by p the perturbation 
pressure nondimensionalized by — P„v^. Then the perturbation velocity potential sat- 
isfies the unsteady wave equation 

V2(j,-M2(^+ ik)2(j) = 0 (2.2) 

9x ^ 

where M=v^/c^ is the free stream Mach number (for simplicity we omit the siibscript 
on M) . The pressure is related to the velocity potential by 

p = “2(-^ + ik) (f). (2.3) 

When the boundary conditions at the upper and lower walls are the same, both 
the pressure and velocity potential can be shown to be antisymmetric functions of 
z ; i . e . , 


p(x,-z) = -p(x,z), (J)(x,-z) = -(|)(x,z). (2.4) 

Since the pressure is assumed to be continuous everywhere in the flow field 
except across the airfoil where it is permitted to be discontinuous (vortex surface) 
and since the pressure is assumed to be continuous at the trailing edge (Kutta con- 
dition) , we may write 
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p(x,0) = 0, |x|>l, (2.5i) 

p(x,0) = - jAp(x), |x|<l, (2.52) 

p(x,0) =0, X = 1, ( 2 . 53 ) 

where Ap is the lifting pressure jump across the airfoil. 

At the surface of the airfoil, the flow must be tangent to the surface. Thus, 

■^1 = w(x) , |x|<l. (2.6) 

z=0 

The purely kinematical boundary condition (2.6) is exactly equivalent to the asser- 
tion that zero mass can flow through the surface of the airfoil. The corresponding 
statement for a closed wall wind tunnel is 


ii =“• 

Z=±T\ 


X <“ I 


(2.7i) 


which in view of (2.3), is equivalent to 


ii , 

Z=±T1 


X <co. 


( 2 . 72 ) 


An open j et boundary is defined mathematically as one on which the perturbation 
pressure is zero (see, e.g., [l,p.47] or [2, p.41]). 


P 4. = 0' 

' Z = ±T) 


X <” . 


(2.81) 


An alternate formulation which is equivalent to (2.8j), 


<l>l2.=+Ti = |xi<”, (2.83) 
is obtained by integrating (2.3). 

We point out that an open jet boiindary differs from free air because in the 
latter case the pressure perturbation due. to the presence of a model will vanish 
in general only at infinity. 
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It has been known since the time of Prandtl [11] (see also Glauert [4, p.43] , Garner 
et al II, p.5] and Kraft [12,p.l]) that the closed wall and open jet boundary con- 
ditions (2.1) and ( 2 . 8 ) cause opposite interference effects in the sense that closed 
walls produce more lift than free air and open jets produce less. The desire to 
exploit these opposing effects to enable an interference free tunnel and to avoid 
the problem of choking has led to the construction of ventilated wind tunnels at 
various facilities. 

The physical flow at a ventilated wind tunnel is quite complex and depends, 
among other things, upon viscous and boundary layer phenomena as well as the detail 
construction of the walls themselves. However, at some distance away from the wall, 
the localized effects of individual slots or holes by which the tunnel is ventilated 
will blend into a more homogeneous effect, thereby permitting the introduction of 
homogeneous boundary conditions . 

Two types of homogeneous boundary conditions have been proposed, the so-called 
slotted wall boundary condition and the porous wall boundary condition. The slotted 
wall boundary is physically based on an accelerative or mass effect whereas the por- 
ous wall boundary condition is based on a viscous effect. However, this simple dis- 
tinction seems not to be generally made. Therefore, to emphasize the underlying 
principles of mechanics involved, we shall adopt in the sequel the more descriptive 
terms mass effect and viscous effect . 

The mass effect boundary condition may be developed heuristically as follows. 
Recall that Euler's vector equation for conservation of linear momentum in a perfect 
fluid is 

(physical dimensions) - ^p = pa. ( 2 . 9 j) 

If one applies this law across a ventilated wind tunnel wall over a distance Z such 
that local effects at the wall may be neglected, there results 

(physical dimensions) ~ tunnel _ ^ plen^ _ p._ ^^ ^ ( 2 . 82 ) 

where the pressure gradient has been approximated by a finite difference quotient 
and where v^^ is the average velocity component normal to the wall. Neglecting the 
difference between the plenum chamber pressure and the free stream atmospheric pres- 
sure (this assumption can be removed [13]), and nondimensionalizing all quantities 
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gives 




z = +n. 


(2.93) 


where 


y = Z/b. 

Regarding the length ratio y as an empirically determined quantity in order to 
account for only part of the wall being closed, we may call y the mass effect 
ventilation coefficient and we call { 2 . 93 ) the mass effect boundary condition . 

An alternate form of ( 2 . 93 ) that is often used, especially in steady flow, is 

!|i+y^ =0, z = +^n- (2.94) 


The mass effect boundary condition ( 2 . 93 ) was apparently first proposed by 
Davis and Moore [14] in 1953. Although their basic analysis does not seem to preclude 
application to other forms of wall ventilation, they were primarily interested in 
slotted wall tunnels and called ( 2 . 93 ) the slotted wall boundary condition. 

Formulas for estimating the mass effect ventilation coefficient y for slotted 
wall tunnels in terms of the slot geometry are given by Davis and Moore [14], 
Guderley [15], Baldwin, Turner and Knechtel [16], Chen and Mears [17] (see [13]), 
Goethert [2 ] , Garner et al [1 ] , Pindzola and Lo [18] , Barnwell [13] and others. 

Recently, Barnwell [19] has concluded that better agreement with experiment 
using the boundary condition (2.9) is obtained if y is determined by correlation 
with experiment rather than if it is calculated using existing theories based on 
slot geometry, etc. 

The viscous effect boundary condition is based on the assumption that the 
wall is porous in the sense that the mass flow rate through the wall is proportional 
to the difference in pressure between the inside of the wind tunnel and the plenum 
chamber; i. e. , 


(physical dimensions) 


P,. -.“P n = pv V , 

tunnel plenum p n 


( 2 .IO 1 ) 


where v is a physical constant having dimensions of velocity and which depends upon 
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the porosity of the wall and the properties of the fluid in the wind tunnel (usu^ 
ally air) . 

If Vp is nondimensionalized by the free stream velocity according to 

_ _ ^tunnel ^plenum , ( 2 .IO 2 ) 

“ V ~ PVr|V 

00 II 00 

and if the difference between the plenum chamber pressure and the free stream atmos- 
pheric pressure is again neglected, then it can be shown that the viscous effect 
boundary condition in the form 


p = +2v^, z = +n (2.IO3) 

— dZ — 

follows. The dimensionless quantity v is called the viscous effect ventilation co- 
efficient . A convenient alternate form of (2.IO3) is 

(-;^+ik)p+v^ = 0, z = +13. (2.IO4) 

oX — OZ — 

The boundary condition described by (2.10) clearly represents a viscous mechan- 
ism in the sense that force in the form of pressure is proportional to velocity. 
Sometimes this is loosely referred to as Darcy's law [20] a la soil mechanics but 
that does not seem to be an appropriate attribution in the present context. Such 
a viscous effect boundary condition for wind tunnels was proposed apparently first 
by Goodman [ 21 ] in 1950 for steady flow. Further investigations along this line 
have been made by Baldwin, Turner and Knechtel [16], Woods [22] and [23] , Drake [24] 
and [ 25 ] > Parkinson and Lim [26] j Ebihara [27], Mokry [28], Kraft and Lo [29] and 
others. Woods and Drake considered walls with finite length porous sections and 
Drake considered unsteady flow but gave no numerical results for subsonic flow. 

In connection with Barnwell's conclusions cited above concerning the difficulty 
of obtaining good agreement with experimental data using the so-called slotted wall 
(mass effect) boundary condition, we note that experimental evidence exists [26] 
which indicates that the viscous effect boundary condition is more realistic for 
slotted wall tunnels with narrow slots (see also [28,p.48]). 

A complete understanding of the mass effect and viscous effect boundary condi- 
tions has yet to be achieved. On the theoretical side, this is partly because the 
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kernel of the integral equation relating downwash and pressure has not been explicit- 
ly calculated for the latter. However, a computationally tractable form of the ker- 
nel for the Tinsteady mass effect boundary condition was presented by Bland [8] in 
1968. This is discussed briefly in section 3 below and further information may be 
found in [5] and [9]. Extensive calculations using Bland's kernel were first 
given in 1978 by us [5] and a rigorous existence and convergence theory for the nu- 
merical solution method was established. While the bulk of our present analysis 
centers around Bland's kernel, we present in section 4 an analysis of the kernel 
for the viscous effect boundary condition for purely porous wind tunnel walls. 

In section 13 we present calculations for steady airloads using the viscous effect 
boundary condition over the full range of the Mach number, height to chord ratio and 
ventilation coefficient v. 

Most of the existing wind tunnel interference theory is founded on the incom- 
plete point of view that the effect of the boundary conditions at the wall can be 
determined by knowing the far flowfield characteristics of the airfoil. While this 
viewpoint has produced useful and simple engineering approximations to angle of at- 
tack and Mach number "corrections", it suffers by neglecting the truly coupled 
nature of interference between the walls and the model. Thus, any rational theory 
of wind tunnel interference must be based on an appropriate boundary value problem 
in which this coupling is explicitly present. The analysis used in this report is 
based on an integral equation method which correctly accounts for such coupling. 

Whether the governing boundary value problem is solved directly via partial 
differential equations or by reduction to an integral equation is just a matter of 
computational method. However, whereas integral equations are more difficult to 
formulate, they enjoy the advantage that the dimension of the space in which the 
problem must be solved is reduced by one. This undoubtedly accounts for integral 
equations being a frequent method of choice in subsonic and supersonic flow, and 
there is now evidence [29] to suggest that in transonic flow as well, integral 
equation methods may be an order of magnitude faster than PDE methods. 
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§3. Bland's integral equation 


Based on the preceding discussion, the boundary value problem for flow in 
a ventilated wind tunnel with the mass effect boundary condition may be completely 
formulated on the lower {or upper) half of the infinite strip 


n = { (x,z) I 

using equations (2.1)-(2.3), (2,5), 
transform pairs 


I X I <“ & 0£z£n } 

(2.6) and (2.9^). 


(3.1) 

Introducing streamwise Fourier 


f (s,z) 


f (x,z) 


1 

/2n 



• oo 

-isx, , . , 

e f(x,z)dx, 

— CO 

. 

ISXJ, , , 

e f(s,z)ds, 

— GO 


(3.2i) 


(3.22) 


with the property that 


■ a 


it follows that 


0 - 3 ^(-^Hst^)p = o, o<z<n, 


l-M 1+M 


p =-j-~ =j e z = 0, 


M ^ + P = 0 , z = Til 

QZ 


(3.23) 


(3.3i) 


(3.32) 


(3.33) 


which can be solved, giving 


p(s,z) 


1 gncr cosh g(n~z)CT + sinh B(n-z)g 
2 /^ gya cosh gria + sinh gna 


^e“"-®^Ap(?)d5, 

-1 


( 3 . 4 i) 
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where 


/ ^ Mk ^ ^ Mk . 

a(s) = As-— )(s+— ). 

Equation (2.6) for the downwash may be restated as 


= /27 

and from equation (2.3) it follows that 


isx,3(j). , 

e (-^) „ds, 

3y y=0 


(3.42) 


(3.5) 


2(s+k) 


(3.6) 


Combining equations (3.4i), (3.5) and (3.6) results in 

.1 


w(x) = 


-1 


BAp (g) is(x-g) io (s) Byg sinh BnP + cosh Bnc 
4 ^ 2(s+k) Byg cosh Bng + sinh Bho 


(3.7) 


Equation (3.7) may be interpreted as an integral equation of the form 


w (x) = — 
4T7 


K(x-g) Ap(g)dg, 


(3.8) 


where the kernel is given by 

isx ig(s) Byg sinh BhC + cosh Bng 


K(x,k,M, Tiry) = 


2(s+k) Byg cosh Bng + sinh Bng 


ds . 


(3.9) 


The inverse Fourier transform (3.9) was evaluated by Bland [ 8 1 following an 
extensive analysis. He showed that 


1 ik , I I . TT . , ,,1+yk tanh kg / ^ 

K(x,k,M,n,y) log lx|+ — (l+sgn(x)) exp(-ikx) 

^ ^ y+i tanh kg 


^ rcsch ^ 4 (exn 

2Bn ^ 2Bn ^x ^ ^ 


2Br/ 


§[log(^ tanh ^) - (exp ^ 


- 1) log tanh , 


(3.10) 
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where 


S{ 6 ) 


^ exp(-An6) 



exp[- 


■ (n- — )it5] } 


(3. Ill) 


s 

n 


[1 


1+(yA„) 


r] [1 + 




( 3 . 112 ) 


= 


= ex 


Mkn, 


(3.II3) 


tan Ajj + yAfj = 0 , 


(3.II4) 


Y 



(3.II5) 


Equation (3.8) is recognized as a Fredholm integral equation of the first kind. 
Its kernel, given by (3.10), identifies it further as a Cauchy singular equation be- 
cause of the dominant — singularity. The remainder of the kernel consists of a 
weaker, integrable logarithmic singularity followed by a continuous part. 

It is well-known (see, e.g. [5,p.ll]) that solutions of such Cauchy singular 
integral equations are not unique unless the auxiliary Kutta condition, given in 
strong form by ( 2 . 83 ) or in weak form by 


lim|Ap(x)|< °°, (3.12i) 

X ^1" 

is imposed. Then it follows also that solutions for Ap possess inverse square root 
singularities at the leading edge 

lim I /l+x Ap(x) I <“ ( 3 . 122 ) 

X -S--1+ 

and square root singularities 


lim I 

x-^ 1" 


Ap(x ) 

y i-x 


<00 


(3.123) 


at the trailing edge. 
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In the case of steady flow in free air, Bland's kernel (3.10) reduces to 
the classical Prandtl-Glauert vortex kernel 

K{x,0,M,«, . ) = ^ (3.13) 

tl 

which consists of the Cauchy singularity alone. In this case, the Sohngen inversion 
formula [30] (see also [31], [32] and [5, §§2,9]) holds, and (3.8) has the closed 
form solution 


(k-o, „-) -It / iSf »<«>«• 

-1 

In view of the above considerations concerning uniqueness, leading and trailing 
edge singularities and the Sohngen inversion formula. Bland changed the unknown in 
(3.8) from the pressure jump Ap to the pressure factor 4* according to 


Ap(x) 


1 /S 

B / 1+x 


iMx) . 


(3.15) 


Then (3.8) becomes 


w(x) 




K(x-?)4^(C)d?, 


(3.16) 


where K is still given by (3.10). The theoretical advantage of (3.16) over (3.8) 
is that (3.16) has a unique solution with the correct leading and trailing edge 
singularities. Without the auxiliary Kutta condition, (3.8) does not have a unique 
solution and is thus not well posed [ 33 , Ch. m , §6. 2] . Since this reformulation of the 
airfoil integral equation into a well posed problem was first made by S.R. Bland 
[ 8 ] , [ 9 ] , we shall refer to (3. 16) as Bland's integral equation. 
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§4. Extension to porous wall wind tunnels 

In a porous wall wind tunnel, the predominant mechanism by which ventilation 
is resisted is viscosity. The governing boundary value problem is given by equa- 
tions (2.1)-{2.6) together with the viscous effect boundary condition ( 2 .IO 4 ). 
Taking streamwise Fourier transforms as in section 3 gives the two point boundary 
value problem 


dz2 


= 0, 0^z_^T|f 


(4.1i) 


P 



-is? 


Ap(C)d?, z 


11 


0 , 


( 4 .I 2 ) 


ii 


V + i(s+k)p =0, z = r|. 
dz 


( 4 .I 3 ) 


Solving (4.1) gives 

- , > 1 Bva cosh B(n-z)c + i(s+k) sinh B(n-z)a -is?, 

= “ 17 ^ Bva cosh Bna + i(s+k) sinh Bna Ap(?)d?, (4.2) 

and relating p to the downwash using (3.5) and (3.6) results in 


w(x) = - 

TT 


^ BAp(? ) is(x-^) icr(s) Bvcr sinh Btcr + i(s+k) cosh Bho 
4 ® 2(s+k) B ^)0 cosh Bncr + i(s+k) sinh Bhcr 

-1 


Equation (4.3) may be interpreted as an extension of Bland's integral equation 
(3.16) to the viscous effect boundary condition ( 2 .IO 4 ). The kernel is given by 


K(x,M,k,n,v) = 


io(s) Bvc sinh Bncr + i(s+k) cosh Bn(7 


2 (s+k) Bva cosh Bto + i(s+k) sinh Bhcr 


ds. 


(4.4) 


and depends upon Mach number, reduced frequency, height to chord ratio and the 
viscous effect ventilation coefficient. 

Technically, both (3.9) and (4.4) are inverse Fourier transforms in the dis- 
tributional sense [34], For this reason, their numerical computation is not straight- 
forward. However, for steady flow the kernel in (4.4) can be obtained easily in 
closed form. 
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Setting k = 0 in (4.4) gives 


K(x,M,n/V) = — 


^isx 3v sinh gns + i cosh 3ns 

3v cosh 3ns + i sinh 3ns ' 


(4.5) 


and introducing a compressible leakage angle C defined by 


tan C = — ~ 
^ 3v 


(4.7) 


results in 


K(x,M,n,v) = 


isx 

e tanh(3ns+i?)ds. 


(4.7) 


For infinite height to chord ratio, n = “, and we observe that (4.7) reduces to 


K(x) = - 


e sgn(s)ds = — , 

X 


(4.8) 


which is the known classical result for steady compressible flow in a free atmosphere. 

The integrals in (4.5), (4.7) and (4.8) do not exist as ordinary Lebesgue (nor 

Riemann) integrals; instead they must be regarded as distributions. However, if we 
write 

K(x,M,ri,v) = — + AK(x,M,n,v) , (4.9) 

X 

then the inverse Fourier transform. 


AK(x,M, n, v) 


_i 

2 


■ oo 

e (tanh 

— OO 


(3vs+ii;) -sgn (s) ) ds , 


(4.10) 


for the incremental part AK of the kernel exists as the Cauchy principal value of an 
ordinary Legesgue integral. Splitting the integral (4.10) into two parts, we obtain 

-e 


AK(x,M,n,v) = i lim( 
e^O 


2 (3ns+i?) 


1+e 


2 (3ns+iC) 


ds 


^-2 (3ns+i?) 


1+e 


-2 (3ns+ii;) 


7 --. ds) , 


and expanding the denominators into geometric series gives 

-e 


AK(x,M,r|,v) = i lim 
e->^0 J 


V , T.n+l (2n3ri+ix) s+2n^i 
I (-1) e ds 


n=l 


i lim 
£->■0 J 


OO oo 


I (- 1 ) 

e n=l 


n+1 (-2n3Ti+ix) s-2nCij 

e ds . 
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For s 7 ^ 0, these series are bounded by an absolutely integrable function, so by 
the dominated convergence theorem, integration and summation can be interchanged. 
Upon taking the limit as e 0, 


■» , 2n?i -2n?i 

. / ^ IV/,, n-1 ,e e 

4K(X,M,„,V) . — I (-1) 

'--I "■^2^ 

which may be written as the imaginary part of a complex series 


(4. Ill) 


1 " 

AK(x,M,ri/V) = — Im ^ ^ — , 


n=l n- 


2 Bn 


( 4 . 1 I 2 ) 


or as a real series 


T ” n sin 2n^-— — cos 2n? 

AK(x,M,n,v) = I (-if—- 

pn 


n=l 




( 4 . 1 I 3 ) 


The last series is the sum of two series which may be summed in closed form [35, 
Nos. 561, 562]. Thus, we obtain the incremental part of the kernel in closed form 
as 


exp 

AK(x,M,n,v) - — ). 

Bn 2 , . -rrx X 

— sinh 

7T 2Bn 


(4.12) 


Adding the Cauchy singularity gives the complete kernel. 


K(x,M,n,v) = 


,Cx. 

^ 4i4h If 
TT 2Bn 


(4.13) 


We point out that, physically, AK represents the interference effect due to 
the presence of the wind tunnel walls since 

lim AK(x,M,n,V) = 0. 

n'>oo 

Thus, we are justified in referring to AK as the interference kernel for steady 
flow . Also, we observe that for a closed wall, (4.13) reduces to 

TT 7TX 

K(x,M,n,“) = csch ^ 75 — , (4.14) 

2Bn 2Bn 

which is a result originally given by Runyan, Woolston and Rainey [36, Appendix] 

(see also Bland [ 9 , p. 839]). 
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The parametric behavior of the interference kernel and the complete kernel 
can be displayed explicitly by writing 

Br|AK(x,M,ri,v) = F(^,t), (4.15) 

BnK(x,M,n,v) =— +f(^,t), (4.16) 

T 


where 


is a function of 


F(?,t) 
two parameters ? 


exp (gT) _ 1 

2 . , IT T 

— sinh — T 
IT 2 

and T given by 


(4.17) 




tan ^ 


1 

Bv' 


T 


X 

Bn* 


Thus, our analysis reduces the original four parameters 


to only three 


x,M,n,v 



Bv . 


(4.18) 


The parameters 

Bn, Bv 

could be expected from the Lorentz compressibility contraction (as pointed out by 
Kiissner for unsteady flow (and by Glauert and Prandtl for steady flow) ) . The para- 
meter 

X 

Bn 

could be expected from the previous result (4.14) for closed wall tunnels and rep- 
resents the ratio of the chordwise influence distance to the contracted height to 
chord ratio. These results were originally derived by Ebihara [27] using the method 
of images in conjunction with compressibility corrections to steady incompressible 
flow (see also Mokry [28] ) . 

The ease of obtaining the interference kernel in terms of elementary functions 
can be attributed to the relative simplicity of steady flow. For unsteady flow, it 
seems unlikely that the inverse Fourier transform (4.4) would admit such a pat re- 
sult. If the singularities were removed, however, the interference kernel for un- 
steady flow in a wind tunnel with porous walls might be expressible as an ordinary 
integral. This integral, once it were obtained, could be attacked via infinite 
series or via numerical integration. For analytic functions, the Laguerre-Gaussian 
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quadrature rule [ 37] is convergent in the sense that weights and nodes are 
tabulated [38] such that 



On the other hand, the uniformly convergent series presented in (4.11) are not 
rapidly convergent. We have found that three decimal accuracy requires tens of 
thousands of terms; such behavior, even with convergence acceleration, is quite 
pronounced in the numerical computation of Bland's kernel (3. 10) - (3. 11) . This 
poses the interesting question: from the standpoint of computational efficiency, 

is it better to calculate the continuous part of the kernel as an infinite series 
as has been customary in the past, or is it better to compute it as a Fourier trans- 
form? 

For future reference, we present tabulated values of the interference kernel 
and of the complete kernel in Tables 1 and 2. Two features may be noted. First, 
it is easy to show from (4.17) that 


lim P(5 ,t) = C- (4.19) 

T^d 


Second, we observe that 

lim F(^,t) = 0 if 0<5<^, (4. 20]) 

lim F(^,t) = 0, ( 4 . 2 O 2 ) 

lim F(^,t) = IT. ( 4 . 2 O 3 ) 

Figures 2 and 3 depict the parametric behavior of the interference kernel and 
of the complete kernel for variable porosity tunnels. In Figure 2, it is seen that 
for ^ = 0 (the closed wall case) , the interference kernel is an odd function of 

streamwise distance. Thus, upstream and downstream distances produce opposing in- 
terference effects. However, the interference downwash, given by 


1 

TT 


1 


-1 



AK(x-e,M,n,v)(J;(?)dC, 


is skewed upstream by the singularity factor 


/ 


1-K 

l+C 


so that upstream pressures have a greater effect on downwash than downstream pressures. 
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Also the interference kernel modifies the complete kernel so that the same pressure 
produces less downwash than in the free atmosphere. Reversing the argrament, the 
same downwash in a closed wind tunnel corresponds to more pressure than in free 
air. As ventilation occurs, the leakage angle 






increases. This skews the interference kernel to the right as shown in Figure 2 
and increases the magnitude of the complete kernel as shown in Figure 3. Thus, 
increasing the ventilation causes upstream pressures to have a relatively greater 
interference effect than downstream pressures, and causes the same pressure to 
produce more downwash than in a closed tunnel. Again reversing the argument the 
kernels shown in Figures 2 and 3 indicate that as the ventilation increases, the 
airloads will decrease for the same downwash. 
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Table 

1. Steadv 

interference kernel 

6 nAK for 

variable 

Dorositv wind tunnels 

X 

T= — 

3n 

5=0® 

C=15° 

?=30° 

5=450 

e=60° 

?=75° 

C=90° 

- 8.00 

.124989 

.124999 

.125000 

.125000 

.125000 

.125000 

.125000 

-7.00 

.142804 

.142849 

.142856 

.142857 

.142857 

.142857 

.142857 

- 6.00 

.166413 

.166614 

.166656 

.166664 

.166666 

.166667 

.166667 

-5.00 

.198780 

.199671 

.199911 

.199976 

.199994 

.199998 

.200000 

-4.00 

.244133 

.247941 

.249278 

.249746 

.249911 

.249969 

.249989 

-3.00 

.305109 

.320465 

.327466 

.330658 

.332114 

.332777 

.333080 

-2.50 

.338078 

.367819 

.383275 

.391308 

.395483 

.397652 

.398780 

- 2.00 

.363985 

.419427 

.452270 

.471725 

.483251 

.490078 

.494122 

-1.50 

.366207 

.463787 

.529767 

.574166 

.604207 

.624492 

.638189 

- 1.00 

.317431 

.474650 

.595656 

.688791 

.760473 

.815644 

.858108 

-.50 

.191725 

.413590 

.608233 

.778995 

.928805 

1.06023 

1.17554 

0 

0 

.261799 

.523599 

.785398 

1.04720 

1.30900 

1.57080 

.50 

-.191725 

.061169 

.349431 

.678006 

1.05254 

1.47944 

1.96605 

1.00 

-.317431 

-.113161 

.152240 

.497066 

.945087 

1.52719 

2. 28349 

1.50 

-.366207 

-.221694 

-.007674 

.309284 

.778689 

1.47387 

2.50340 

2.00 

-.363985 

-.270395 

-.112406 

.154295 

.604510 

1.36451 

2.64747 

2.50 

-.338078 

-.280850 

-. 170735 

.041147 

.448845 

1.23333 

2.74281 

3.00 

-.305109 

-.271430 

-.197562 

-.035548 

.319792 

1.09915 

2.80851 

4.00 

-.244133 

-.233282 

-. 202359 

-.114239 

.136872 

.852451 

2.89160 

5.00 

-.198780 

-.195485 

-.183282 

-.138102 

.029176 

.648516 

2.94159 

6.00 

-.166413 

-.165447 

-.160800 

-.138445 

-.030906 

.486406 

2.97493 

7.00 

-.142804 

-.142528 

-.140798 

-.129990 

-.062435 

.359791 

2.99874 

8.00 

-.124989 

-.124911 

-.124278 

-.119133 

-^07735^ 

_,_2.61_87_1 _ 

3.01659 


Table 2. Steady complete kernel BnK for variable porosity wind tunnels 


H 

II 

II 

0 

0 

?=15° 

C=30° 

5=450 

C=60° 

C=750 

?=90° 

- 8.00 

-.000011 

-.000001 

-.000000 

-. 000000 

-.000000 

-.000000 

-.000000 

- 7.00 

-.000053 

-.000008 

-.000001 

-.000000 

-.000000 

-.000000 

-.000000 

-6. 00 

-.000254 

-.000053 

-.000011 

-.000002 

-.000000 

-.000000 

-.000000 

-5. 00 

-.001220 

-.000329 

-.000089 

-. 000024 

-.000006 

-.000002 

-.000000 

-4.00 

-.005867 

-.002059 

-.000722 

-.000254 

-.000089 

-.000031 

-.000011 

-3.00 

-.028224 

-.012868 

-.005867 

-.002675 

-.001220 

-.000556 

-.000254 

-2.50 

-. 061922 

-.032181 

-.016725 

-.008692 

-.004517 

-.002348 

-.001220 

-2.00 

-.136015 

-.080573 

-.047730 

-.028275 

-.016749 

-.009922 

-.005878 

-1.50 

-.300460 

-.202880 

-.136991 

-.092501 

-.062459 

-.042175 

-.028478 

-1.00 

-.682569 

-.525350 

-.404344 

-.311209 

-.239527 

-.184356 

-.141892 

-.50 

-1.80828 

-1.58641 

-1.39177 

-1.22101 

-1.07120 

-.939766 

-.824462 

0 

+00 

+ C0 

+ 00 

+00 

+ 00 

+00 

+00 

.50 

1.80828 

2.06117 

2.34943 

2.67801 

3.05254 

3.05254 

3.96605 

1.00 

.682569 

. 886839 

1. 15224 

1.49707 

1.94509 

2.52719 

3.28349 

1.50 

. 300460 

. 444973 

.658993 

.975950 

1.44536 

2. 14053 

3.17007 

2.00 

.136015 

. 229605 

.387594 

.654295 

1.10451 

1.86451 

3.14747 

2. 50 

.061922 

.119150 

.229265 

■ .441147 

.848845 

1.63333 

3.14281 

3. 00 

.028224 

.061903 

.135771 

.297785 

.653126 

1.43249 

3.14185 

4.00 

. 005867 

. 016718 

.047641 

. 135761 

.386872 

1.10245 

3.14160 

5.00 

.001220 

.004515 

.016718 

.061898 

.229176 

.848516 

3.14159 

6.00 

.000254 

.001220 

.005867 

.028222 

.135761 

.653073 

3.14159 

7. 00 

.000053 

.000329 

.002059 

.012867 

.080422 

.502648 

3.14159 

8.00 

. 000011 

. 000089 

.000722 

.005867 

.047641 

.386871 

3.14159 
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Steady interference kernel for variable porosity wind tunnels 



Figure 3. Steady co mplete kerne^ f or vari a ble porosity wind tunnels 
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§5. Fourier theory of airfoil polynomials 


In our previous work [ 5 ] , two sets of orthogonal polynomials were used in 
the solution of Bland's integral equation. These polynomials are computationally 
well suited for unsteady flow in ventilated wind tunnels, probably because the air- 
loads for subcritical Mach numbers and nonresonant frequencies are well-behaved 
parametric extensions of the steady free air case which they solve so elegantly 
in closed form. From a theoretical viewpoint, they have unique mathematical pro- 
perties which provide efficient use of the abstract geometrical structure of the 
underlying Hilbert spaces of downwash and pressure functions. This section briefly 
siammarizes some of their features. For additional information, see Bland [8 ] and 
[9 ], ivanoff [7, Ch. 2 ] , or Fromme and Golberg [5 ]. 

Consider the case of steady flow in free air. A completely symmetrical theory 
between downwash and pressure will now be presented for this special case. To do 
this, (3.14) and (3.16) are viewed as a solution pair 


-1 

1 

1 


'izi JL 

1+? X-C 

/TTe 1 


IMC) dC = W(x) , 


w(?) dC = li-(x) . 


1-C ?-X 

This may be written more briefly in operator notation as 


where 


a (x) 
n 


Y (x) 
n 


^ cos ((n-^) cos ^x) 
cos(-^os“lx) 

^ sin ( (n-y) cos ^x) 
sin(^os“^x) 


(5.1i) 


( 5 .I 2 ) 



H 41 

= W, 

(5. 

•2l) 

h' 

■1 w 

= 

(5, 

. 22 ) 

independent 

solution pairs to (5.2) is given by 



H Y = 

a , 

n = 

(5. 

. 3i) 

n 

n 



a = 
n 


n = 1, . . . ,®, 

(5. 

.32) 


(5.4j) 


(5.42) 
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These functions satisfy the following recursion formulas 


a^(x) = a2W = ^(-l+2x) , . . . , 

“n+2^^> 

= 2x a (x) -ct (x) , 

n+1 n 

(5.5i) 

^ = ^(l+2x),.... 


= 2x Y . T (x)-Y (x) , 
n+i n 

(5.52) 


and are therefore polynomials. We call a and y the nth downwash and pressure poly- 

n n — 

nomials , respectively. Collectively they are called airfoil polynomials . Their n-1 
zeros are given by 


a^Cx^) 

= 0; 

n 

X . = 

1 

2iir 

2n-l' 

i = 1,. 

• • r ri“ X t 

n>2 , 

( 5 . 61) 

Y (C") 
n 1 

= 0; 


2iTT 

cos ? 

2n-l 

i = 1,. 

- • f n-1; 

n>2. 

(5.62) 


and are interdigitated according to 


-1 < 


’n -1 


n 

< < 




X < 1. 

n -1 


(5.63) 


One of the more striking properties of the airfoil polynomials is that they are 
orthogonal with respect to reciprocal weight functions with leading and trailing edge 
singularities : 


a. (x)a (x)dx = 
/ 1 -x m n 


5 , 


J ^ / 1+x 


Y (x)y (x)dx = 6 


(5.7j) 


( 5 .? 2 ) 


This leads to two generalized Fourier series representations, one for downwash, the 
other for pressure. To see this, define two complex Hilbert spaces 


= 




{f; [- 1 , 1 ] 

{f: [-1,1] >C| 



jf (x) I dx<«>}, 
2 

I f (x) I dx«”} , 


(5. 81 ) 


(5.82) 


called downwash space and pres sur e spac e , with respective inner products 
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<f,g> = 

a 


<f,g>^ = 


n /i^ 
J / l-x 
-1 

/lEE 

J / l+x 


f (x)g(x)dx. 


f (x)g(x)dx. 


(5.91) 

(5.92) 


and norms 


1 f II = /<f ,f> , 
' " a a 


lf|| = Af,f> , 
Y Y 


(B.lOj) 


( 5 .IO 2 ) 


It can be shown that {a } and {v } are complete orthonorraal bases for I- and L , 

n n a y 

respectively. Therefore, an arbitrary function in either space can be represented by 
generalized Fourier series using airfoil polynomials; i.e.. 


feL^ => 

f = ) <f,a > a , 

mam 

(5. Ill) 

a 


m=l 



00 


=> 

l-tl 

11 

A 

l-h 

V 

( 5 .II 2 ) 

Y 

1 ^ f ^ 

m=i 



Referring to the solution pair given by (5.1), we see that for steady flow in 
free air. 


wsL^ => i^eL^ & i); = T <w,a > y , 
a Y mam 

m=l 

00 

\jjo => w€ & w = J <\jj,y > a 

^ A# ^ ^ ' m *v 1 


m=l 


my 


(5.121) 

( 5 . 122 ) 


Thus, in this fundamental special case, the pressure Fourier coefficients and the 
downwash Fourier coefficients are exactly equal; i.e., 


<{lj,y > = <w,a > , (5.13) 

my n* c( 

It 

thereby providing an elegant and computationally powerful reformulation of the Sohngen 
inversion formula. 

Another way of viewing the result (5.13) is that if t|j and w satisfy (5.1), then 
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they have the same length, or norm, in their respective space; i.e.. 


Hi/) = w => II i/)ll^ = II w 11^ 


(5.14) 


In general, if T is a linear transformation between two Hilbert spaces and H 2 , 


and if we define the operator norm of T as 


(5.15) 


||t|| = sup II Tx| 

11 x 11^=1 

where ||-||. is the norm on H. i=l,2, and if 
1 1 / 


(5.16) 


IMI = 1 - 


(5.17) 


then we say T is an isometry. By (5.13), H is an isometry; i.e.. 


||h|| = 1 , ||h mi = 1- (5.18) 

An isometry is an abstract generalization of a length preserving transformation and 
enjoys generalized properties of orthogonal, or unitary finite dimensional transform- 
ations. If T is defined in general as above, then the adjoint of T is the operator 

T*:H2 Hi (5.19) 


such that 


<T*x,y>i = <x,Ty >2 (5,20) 

for every xeHi and ycH 2 , where <.,.>^ are the inner products on H^, i=l,2. It follows 
(see [5] for details) that the adjoint of H is its inverse; i.e., 

H* = H~^. 


(5.21) 
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Consequently , 

fcL^ & geL^ => <f,Hg> =<H”lf,g> , (5.22) 

a y a Y 

which is useful in reverse flow formulations. 

Certain specific Fourier expansions in terms of airfoil polynomials will be 
useful in the sequel. Let 


log^ ? = log lx-? I . 


Then it has been shown [ 8 ] that 


<log ,a > - 

X n oi 


Y (x) 
n+1 

2n 


(log2-— )Y (x) , n=l, 
2 n 


Y (x) Y (x) Y„ T 
n+1 n , n-1 


<loq ,CX ^ = 

X n a 

2n 

2n(n-l) ' 2 (n-1) ' 

n>z , 

<log ,y > = 

X n Y 

a (x) 

n+1 

+ (log2--^)a (x) , n=l 
2 n 


2n 

1 


ct (x) 

a (x) a T (x) 


<log ,y > = 

^x n Y 

n+1 

n n-1 

n>2 . 

2n 

2n(n-l) 2 (n-1) ' 


(5.23) 


(5.24j) 


(5.242) 


(5.243) 

(5.244) 


These results are used in the solution of Bland's integral equation. For problems in- 
volving flaps, jump functions and their powers given by 


<x> 


k 


(max (0,x) ) 


k 


(5.25) 


are employed. We briefly indicate the method of computing their airfoil polynomial 
Fourier coefficients. To evaluate 


use the changes of 


Since 


k 1 /l-x , ,k 

«x-a> ,Y > - - I /t-T (>^-a) 


variables represented by 

X = cos 9, a = cos6_ 


f l-x 

1+x 


2 sin^ — 


Y (x)dx, 
n 


(5.26) 


sin 0 
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and 


sin (n- 
" sxn - 


0 


it follows from the elementary identity 


2 sin a sin 0 = cos(a-S) - cos(a+B), 


V 

(cos 9-cos 9^) (cos (n-1) 0-cos n0)d0. (5.27) 

The advantage of (5.27) over (5.26) is that the integrand 

F , (0) = (cos 0-COS 0 ) (cos (n-1) 9-cos n9) 
nK a 

is an analytic function everywhere on [0,9^] and therefore (5.29) can be computed ac- 
curately and conveniently using ordinary Legendre-Gaussian quadrature. Thus, 


that 


k 1 

<<x-a> ,y > = 7 = 

n Y 


<<x-a> ,Y 


n y 




I w^- 
1 


F , (x” 9 ) + 
nk 1 a 


-N' 


where and x^ are weights and nodes of the N point Legendre-Gaussian quadrature rule, 
and where 


E 


N 


2N+1 ^ 

2 < N!) (2N) - 

(2N+1) (2N! ) 3 nk ' 


O<0<0 

a 


is the error [37, §8.5]. The integrals (5.26) can be determined in closed form, although 
they are increasingly ciunbersome for larger values of k. The following closed form in- 
tegrations are noted. 

<<x-a>®,y > = ;= (9 -sin9 ), (5.28i) 

1 Y /ir a a 

sin (n- 1) 0 sin n 0 

«x-a>0,Y^> = ^ ^ , n> 2 , ( 8 . 282 ) 

n Y /7f n— J_ n — ^ 

. 0 sin 20 

<<x-a>,Yi> = - 7 = (sin 9 ^ - ^) 

1 Y A a 2 4 ' 


1 

- cos 
/tt 


0 ( 0 -sin 0 ) , 
a a a 


(5.29i) 
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, . sin 9 sin 26 sin 30 

, , . . 1 ,0a a . a a, 

- /; 'T - ^ — 


sin 29 


- ^ cos 9a (sin 9_ 

VTT 


-), 


, sin(n-2)0 sin(n-l)6 sin n0 sin(n+l)0 

i , a a a a, 

«x-a>,Y^> = ^ ) 


'n Y /jf 2(n-2) 


2(n-l) 


2n 


2(n+l) 


(5.292) 


sin(n-l)6 sin n6 

;= cos 6 ( ; — - ~ — — — , n>3. 

/tt a n-1 n — 


( 5 . 293 ) 


We also note for future reference the Jacobi-Gaussian quadrature formulas [37] 
■1 N 

-1 


7 j /S ^ 

n= 
N 

f(x)dx = J 


1+xi"^^ . N+1 

— f (x ) +E 


n=l N+ j 


'N 


(5.30i) 


^ /l=2i _ V 


tl n=l N+- 


N ' 


(5.30^) 


N+1 N+1 

where the nodes and x. are the zeros of the airfoil polynomials a„ , and y 

1 1 N+1 N+1 

as given by (5-6) , and where the error for continuous functions is proportional 

to the 2Nth derivative of the integrand at some point in the interval (-1,1). 

Sometimes it is desirable to convert from one airfoil polynomial basis to the 


and 


to (5.5) 

t 

it 

is 

easy 

to show that 


“1 

= 

Yi' 




(5.31i) 

“2 

= 

2 Y 1 

+ 

Y 2 ' 


( 5 . 3 I 2 ) 

“3 

= 

2ll 

- 

2y2 

+ Y3^ 

( 5 . 3 I 3 ) 

«4 

= 

2Yi 

+ 

2y2 

- 2y3 + Y4' 

( 5 . 3 I 4 ) 

“5 


2Y1 

“ 

2y2 

+ 2y3 - 2y4 + Y5' ■ • • 

( 5 . 3 I 5 ) 

Yl 

= 

“ 1 ' 




(5.32j) 

Y2 

= 

2a-^ 

+ 

0t2 ' 


(5.322) 

Y3 

= 

2 a 1 

+ 

2tt2 

+ 03 . 

( 5 . 323 ) 

T4 

= 

2aj 

+ 

2a2 

+ 203 + 

( 5 . 324 ) 

^5 

= 

2aj 

+ 

2tt2 

+ 2c 3 + 2tt4 + U5 , . . . 

(5. 325 ) 

30) , it 

appears 

that 



OL , = 

n+1 

^n+1 

+ 2 

n 

T, 

m=l 

, - . m+1 , „ 

(-1) Y ; n=l,2, 

m 

(5.33j) 

\+l 

= a 

n+1 

+ 2 

ri 

LI a r n = 1,2,... 
m=l m 

(5.332) 
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and in matrix notation that 


■ «1 


’ 1 0 0 0 0 


. . 

Yl 

«2 


-2 1 0 0 0 


Y 2 

“3 


2-2 1 0 0 


Y3 

«4 

= 

-2 2-2 1 0 


Y4 

“5 


2-2 2-2 1 


Y5 

Yl 


- 

1 0 0 0 0 


ai 

Y2 


2 1 0 0 0 


«2 

Y3 


2 2 1 0 0 


“3 

Y4 

= 

2 2 2 1 0 


«4 

Y5 


2 2 2 2 1 


“5 

. 


. 




or, more briefly 


{a } = [G ] {a }. 
n mn n 


(5.34i) 


( 5 . 342 ) 


(5.35i) 


{y } = [ G ^ Ha }. 


( 5 . 352 ) 


The problem of computing downwash from discrete displacement data will entail 
interpolation procedures and, differentiation as well. In any ascending polynomial 
basis, transformation matrices such as (5.34) will be lower triangular and when 
differentiation is performed they will be lower subtriangular. By differentiating the 
polynomial expressions (5.5) and recombining, one easily obtains 
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= 0, (5.36i) 

= 2a^, ( 5 . 362 ) 

03 = 2ai + 402 , (5.363) 

04 = 4ai + 2 o2 + 603, (5.364) 

05 = 4o]^ + 6 o2 + 203 + 804,... (5.365) 


y1 = 0, 

Y2 = 2Yi, 

Y3 =-2Yi + 4y2, 

Y4 = 4yi - 2y2 + 6y3' 

Ys =-4Y1 + 6 y2 - 2y3 + 8y4/--- 


From this, it appears that in general 


n+1 n 


n 





a ; 

n = 

1 , 2 , 

m=l 

m 



n 

(-1) 

n+m 


^1 

m=5l 

m 


and 






" 

ai 


0 0 0 0 0 . . . 



02 


1 0 0 0 0 . . . 


“2 

«3 


1 2 0 0 0 . . . 


«3 

04 


2 13 0 0 


Ct 4 


= 2 






2 3 14 0 


«5 

. 


. 


- 




r *1 


yz 


Y 3 


^4 


^5 


. ■ . 


(5.37i) 

(5.3?2) 

( 5 . 373 ) 

( 5 . 374 ) 

( 5 . 375 ) 


(5.38^) 

(5.38^) 


(5.39i) 


(5.392) 
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Thus we may write the infinite matrix equations 

{a }' = [d“ ]{cx }, 
n mn n 


(5.40i) 


{y }' = [D^ ]{y 

n mn n 


( 5 . 402 ) 


where the differentiating matrices [D* ] and [d"'^ ] are as shown in (5.39). For higher 

mn mn 

derivatives, we have 


{a } 

n 


(k) 


[d“ ] }, 

mn n 


(5.41i) 


{y = [d'^ ] }, 

n mn n 


( 5 . 412 ) 


where (k) refers to the order of the derivative on the left hand side and to the power 
of the differentiating matrix on the right hand side. 

Formulas expressing derivatives of the airfoil polynomials in terms the airfoil 
polynomials can also be derived in general by making use of their Fourier series pro- 
perties. ^ Since the derivative of a polynomial is a polynomial of one degree less, we 
have • n-1 

a' (x) = y <a',ct > a (x) , 
n n m a m 

m=l 

where the Fourier coefficients are given by 


y,' ,a > = [ ex' (x)a (x)dx. 

n m a J / 1-x n m 


Referring to equation (5.4), the substitution 


X = cos 6 , 

together with the use of elementary trigonometric identities leads to 

IT 

1 1 I , , 1, sin (n+m-1) 6+sin (n-m) 0 

<a ,cc > = — ( (n-— ) ; s 

n'ma '''■Jo 2 sin 6 


^ cos (m+n-1) 6+cos (n-m) 6 
2 1+cos 6 


) d6 , l<m<n-l. 


(5.42) 


^S.R. Bland, 


private communication. 


dated 3 August 1977. 
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The tabulated integrals [39,p.366] 


_1 

TT 


sin n6 
sin 9 


de 


0 if n is even. 


1 

TT 


TT 
^ 0 


sin n6 
sin 0 


de 


1 if n is odd. 


1 

TT 


COS n9 
1 + a cos 0 


d6 = 


, /l-a^- 


/l-a^ 




a <1 


permit (5.42) to be calculated in general. Since for every pair (n,m) of integers, 
precisely one of n+m-1, n-m is odd, it follows that 


f TT 

1 sin (n+m-l) e+sin (n-m) 0 

.A clU 

TT J Q sxn 0 


1 . 


By evaluating the limits 


lim 

a->-l 


1 

TT 


COS (m+n-1) 6+cos (n-m) 9 
1+a cos e 


de 


= lim 
a->l“ 


yj^-l) m+n-1 


/z TT m+n-1 

vl-a^ a 


/l^ 


one obtains 


1 

TT 


•IT 

0 


cos (m+n-1) +COS (n-m) 9 
1+cos 9 


de 


(-1)" "*(2m-D . 


In this manner, the following general formulas result. 


<a' ,a > =0 if m>n, 

n m a — 

<«' ,a > = n+m-1 if m<n & n-m is odd, 

n m a 

<a',a > = n-m-2 if m<n & n-m is even, 

n m a 

<Y' , Y > =0 if m>n, 

n m Y ~ 

<y',y > = n+m-1 if m<n & n-m is odd, 

n my 

■^y',Y > = n+m+2 if m<n & n-m is even, 

n m Y 


(5.431) 

( 5 . 432 ) 

(5.433) 

(5.434) 

(5.435) 
(5.43s) 


which verify (5.39). 



Although it appears that Bland was the first to make extensive application of 
the properties of the airfoil polynomials in unsteady flow problems, their use for 
solving singular integral equations goes back at least 40 years. In this regard 
we make the following historical observations. 

First it can be shown that they are suitably renormalized Jacobi polynomials. 
Multhopp [40] (see also [41] ) utilized Jacobi polynomials as early as 1938 to per- 
form chordwise numerical integration on the problem of steady three dimensional flow 
over wings. This is the origin of the successful technique of interdigitated col- 
location and quadrature points used later by Hsu [42] and others for unsteady flow, 
and in the special case of one collocation point it reduces to collocating the down- 
wash at the three-quarter chord. 

The fact that H is unitary was apparently first shown by Akheizer in 1945. We 
became aware of this fact through the recent publication of Ivanov's book [7, p.l33]. 
This result was obtained using the fact that Ha^ = This important property occurs 

as a particular case of a result given by Tricomi in 1951 [31] (see also [43] ) who in 
turn refers to Szego's book, first published in 1939 [44]. 

In the western aerodynamics literature the specific form of the pressure poly- 
nomials first appeared in a 1967 jaaper on unsteady narrow channel flow by Bland, 

Rhyne and Pierce [10] and their properties were further developed in [8] . (See 
also [9] . ) 

In view of the above observations it seems appropriate to name the transforms 

H 

and 

h"^ = H* 

as the Akheizer-Bland transforms. 
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§6. Solution bv collocation of the integral egii.atipji 


Assume WE and il/e . 

a Y 

operator notation as 


Let Bland's integral equation (3.16) be written in 


Tif) = w. 


( 6 . 1 ) 


Since the pressure polynomials {y } are a complete orthonormal basis for the pressure 

2 ^ 
space we have 


'P = I 


n=l 


n Y u 


( 6 . 2 ) 


Because of this, we seek an approximate solution to (6.1) of the form 


'P = J a . Y 

't'u L 

n=l 


(6.3) 


where {a } , are to be determined. Substituting (6.3) into (6.1) gives 

nN n=l 


N 

y a Ty + r = w, 
nN 'n N 

n=l 


(6.4) 


where 


r = w - M = T(tjj-\p ) 
N N N 


(6.5) 


is the resulting downwash residual error. The general collocation equations are 
that 


r (x ) = 0; m = 1 , . . . , N 
N m 


( 6 . 6 ) 


where { x are the collocation points. Equation (6.6) may be written in matrix 

m 1 

form as 


[C ] {a } = {w(x ) }, 
mn nN m 


(6.7) 


where [C ] is the NxN collocation matrix given by 
mn 


C = (Ty ) (x ) . 
mn n m 


( 6 . 8 ) 


To compute the collocation matrix, it is convenient to split the operator T into 
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three parts. Referring to (3.10), we may write 


T = H + K + K , 
Xi c 


(6.9) 


where H is the Akheizer-Bland transform, where K is given by 

Jj 


(K f) (x) = ^ 

L Trg2 


-1 


1+C 


log (x-5 f(5)d?. 


( 6 . 10 ) 


and where is the remaining part of the transformation with continuous kernel. 


(K f) (x) = -- 

C TT 


n 


^ K^(x-5)f(C)a?. 


(6.11) 


In Bland's collocation method, the collocation points Xm are selected as the zeros 
of the downwash polynomial The unitary and logarithmic parts are then 

computed exactly by (5.3) and (5.24), and the continuous part is computed approxi- 
mately by (5.30) using N point Jacobi-Gaussian quadrature. Thus, 


! N+1. ik 

C = a (x )- — ^ <log N+1, 
mn n m 

m 


Y > 
n Y 


N 

+ I 

i=l 


1 

N+j 


^ , N+1 ^N+1, ,^N+1, ,^N 

K (x -5. )y (C. )+E . 

cm 1 n 1 mn 


( 6 . 12 ) 


The collocation solution is then completed upon solving (6.7) using (6.12) and the 
known downwash. 

For closed wall wind tunnels. Bland [9] observed that (6.12) produced rapid 
convergence with N for smooth downwashes but offered no proof of convergence. In 
[5] , we observed similar convergence characteristics for ventilated tunnels and 
established a rigorous proof of convergence of the method based on a three-way 
equivalence between collocation, complex least squares and Galerkin's method. Spe- 
cifically, we showed that solving (6.7) using (6.12) gives 


lim a 
N-x» 


nN 


Vy 


(6.13) 


under very general conditions. In section 10 below we present a direct proof of 

convergence with sharper error estimates and discuss the computational problem of 

weak convergence in the presence of flaps. In section 11, we discuss various 

methods of accelerating convergence and show that a 500 fold reduction in error 

can be achieved for flaps with approximately a 20% increase in computing time. 

N 

In section 12 we demonstrate an improved method of computing for unsteady 

flow which reduces quadrature errors in (6.12) by approximately 1000. 
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§ 7. Computation of airloads 


7.1 Basic equations 

Once the Fourier coefficients known, it is an easy matter to 

compute pressures from the expression 


Ap(x) = — / r— y <ip,y> Y (x) . 
B / 1+x ^ ' Y n 

n=l 


(7.1) 


Let the lift and moment coefficients be defined as in Ashley [45,p.53] according to 


L=— pv^SC, M=— pv^ScC , 

2 oo CXI L 2 °° “ M 


(7.2) 


where L denotes lift, M denotes pitching moment about the quarter chord, positive 
in the direction of increasing angle of attack, S denotes planform surface area, 
and where c denotes mean aerodynamic chord. Then for an airfoil the lift and moment 
coefficients and center of pressure reduce to 




Ap (x) dx 


2i/tt 

T” Vy' 


,1 


~ 4 


-1 


1 "[]■ 

(x+— )Ap(x)dx = - t 


(7.3) 


(7.4) 


’'CP 



(7.5) 


These particularly simple formulas involve only the first two pressure Fourier co- 
efficients because of the orthogonality properties of the airfoil polynomials, and 
are among the most accurately computed quantities in the TWODI program. 

The aerodynamic work matrix is useful in the solution of aeroelastic problems. 
Let the displacement function be expressed formally as 


CO 

h(x,w) = y q^(u)h^(x) (7.6) 

m=l 


where 

sense 


h are displacement basis functions and q^ are generalized coordinates in the 
of Lagrange's equations. Then the components of the aerodynamic work matrix 
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are given by 


or equivalently by 


A 

mn 


g [1 _9h_ 9Ap „ 

i, 


A 

mn 




h (x) ilj (x)dx, 
m n 


(7.7) 


(7.8) 


where \J; represents the pressure factor corresponding to the displacement function 

h . The components of this matrix depend upon the particular deflection basis used 
n 

(which in practice is often selected as the in vacuo vibrational eigenfunctions)^ 

We will assume that all displacement functions and their derivatives belong to 

2 

that is, each h has a downwash function in L Expanding both ij) and h in the 
n ct n n 

pressure Fourier series, and integrating, (7.8) becomes 


[A 

mn 


] 


— [<h ,y >] ^ 1* 

2 m n m n Y j 


(7.10) 


where * denotes the complex conjugate transpose. 
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7,2 Energy properties of the airfoil polynomials 

In our earlier report [5] , we compared the results from the TWODI program 

n 

against the exact Sohngen and Kussner-Schwarz solutions, and made a precise deter- 
mination of the accuracy of TWODI. Originally, we used the downwash polynomials 
as a basis for the deflections and calculated the integrals for the aerodynamic 
work matrix using numerical integration. We have since observed that if one 
selects the pressure polynomials as the displacement basis, then the aerodynamic 
work matrix in steady free atmosphere flow is upper triangular with zero diagonal 
elements. To see this, choose the displacement basis functions as the pressure 
polynomials 

h = y ; n = 1,2,... (7.11) 

n n 


Then the downwash functions 


w 

n 


(7.12) 


are readily found to be 

wi = 0, (7.13i) 

W2 = 2ai, (7.132) 

w = 6 ai , + 4 « 2 » (7.133) 

3 

w = 12ai + 10a2 + 6a3, (7.I34) 

w = 20k2^ + 1802 + 14o 3 + 804,..., (7-13r ) 


using (5.33) 
the above w 


and (5.34). Using (5.12), the pressures Ap^ 
are given by 

‘p' -lyil <»'■ 

An - - (3 Yi-^2y2) . 

AP4 = (6Yi+5Y2+3Y3) . 

= f /l^ (10yi+9y2+7y3+8Y4) , . . . 


corresponding to each of 

(7.14^) 

(7.14^) 

( 7 . 143 ) 

(7.14^) 

(7.14^) 
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Therefore, by inspection, the lift and pitching moment coefficients and center of 
pressure are given by 



0 


L 

0, 




(7.15i) 

^L2 " 

4/ir 
e ' 

^M2 

0, 

X2 = 

.2500, 

(7.152) 


12^" 
6 ' 


4 /t 
6 ' 

^3 = 

.5833, 

( 7 . 153 ) 

^L4 " 

24 /iT 
6 ' 


io/tT 

B ' 

X4 = 

.6667, 

( 7 . 154 ) 


40/t 

B ' 


is/m 

B ' 

X5 = 

.7000, . . . 

(7.15s) 

and the aerodynamic work matrix is 









0 

1 3 

6 10. . 







0 

0 2 

5 9. . 

• 




[A ] ^ 
mn 

= 

0 

0 0 

3 7. . 



(7.16) 




0 

0 0 

0 4 







0 

0 0 

0 0 


1 



Because the components A represent the work done on the structure as it de- 

mn 

forms in mode m against the pressure due to mode n, it follows that for steady free 
atmosphere flow, zero aerodynamic work is done on the structure as it deforms into 
the pure shape of any of the pressure basis functions . The first column is 

trivially zero because the pressure is zero in steady vertical translation. The 
second column A^^ represents the work done as the structure deforms in its various 
modes against the pressure due to the second mode. In the second pressure mode, 
only the first deflection mode produces work on the structure, and since the second 
displacement mode represents a flat plate at uniform angle of attack, the element 
A ^2 corresponds to the lift coefficient derivative The diagonal element 

*22 work done in pitching the airfoil about the node of mode 2; since 

this node and the center of pressure are both at the quarter chord, no aerodynamic 
work is done. Similar interpretations involving the flexible modes may be made for 
columns 3,4, .. . 
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7.3 Transformation properties of tmsteady airloads 

If the unsteady airloads are knovm with respect to one set of structural (dis- 
placement) basis functions, then they can be computed with respect to any other set 
of basis functions, provided the transformation from one basis to the other is known. 
This is a direct consequence of the invariance of airloads under change of basis. 

To see this, let {h } and {h } be two bases for the Hilbert space L? of airfoil 
mm h 

displacements and let their corresponding generalized coordinates. 

Then 

h{x,0)) = ^ q^(w)h^(x) = ^ (7.17) 

m=l 


m=l 


(Throughout this subsection, we will regard all infinite series expressions in a 
formal sense.) Due to the assumed linearity of the problem, the pressure factor de- 
pends linearly upon the displacement and may be expressed as 


I 

m=l 


v» 


00 




(7.18) 


where 

larly. 


il) and ij) represent the pressure factors due to h and h , respectively, 
mm mm 

the lift and pitching moment coefficients are given by 


Simi- 


00 00 


= I ^m 


t gm 


1— 1 

m=l 

CO 

OO 


C — 3 ' 
II 


^ q 

^Mm' 


m=l 


(7.19) 

(7.20) 


where Cj^ and 
due to h^ and h^, 
given by 


and Cj^ and C 
respectively. 


1 ^ represent the lift and pitching moment coefficients 
The components of the aerodynamic work matrix are 


A 

mn 



h dC, 
m n 


A 

mn 



h tp d?. 
m n 


(7.21) 


To establish the transformations between airloads in the {h } basis and the 

m 

{hjjj} basis, we utilize the transformations (assumed known) 


OO 

h = y H h , h 
m mn n m 

n=l 


OO 


I 

n=l 


H h 
mn n 


(7.22) 


between these bases, 
one another 


Clearly these transformation matrices are the inverses of 


m=l m=l 


(7.23) 
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Substituting (7.22) into (7. 18) - (7. 20) and equating coefficients of generalized 
coordinates gives the transformations between pressures and section coefficients, 


= y H Ip f = y H >p , 
tn mn ^n m mn n 

n=l n=l 


(7.24) 


Ct = y H Ct , Ct = y H 
n=l 


T — I 

mn ^n 

n=l 




(7.25) 


'Mm 


= I H 

n=l 


mn 


'Mn ”mn ^M^- 

n=l 


(7.26) 


Combining (7.27) and (7.21), we immediately obtain 


•1 

-1 


00 CD 


^ dC = Z I H 


1+^ m n 


k=l 1=1 


00 00 


mk 


n i_r - ^ 

/ lit *\ « V' 

1 


'izl 

1+C 


\ ^n V J \ V- 

k=l 'C=l 1.1 


(7.27i) 


(7.272) 


Thus, the components of the aerodynamic work matrix transform according to a dOngfUL- 

dnc.0. t/ianA ofmcution : i . e . , 



CO QO 

^ A A 


00 CO 

V C 


A 

mn 

1 1 
k=l 1=1 

V V V' 

A 

mn 

j 

ii Ji 

(7.28) 

transformations may 

be stated in 

matrix 

notation as 




} = [H Hijl } 

f {^ } = 

= [H ]{i)) }, 

(7.29) 


m 

mn n 

n 

mn n 




= tH ]{Ct } 
mn -‘-'n 


- 

(7.30) 



= [H ]{C» } 
mn Mn 


- 

(7.31) 

[A ] 
mn 

= [H ] [A ] [H 
mn mn mn 

[A ] = 
mn 

[H 1 [A ] [H 
mn mn mn 

(7.32) 


The above results are valid for an arbitrary pair {hjjj}, {h^j,} of bases for Hj^. 
The special case = {Ym^ ds of particular interest. Since in this case 


H = <h ,Y > , 
mn m n Y 


(7,33) 
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it follows that 


} - [<h > ] {ijj 

m m n y n 

(7.34) 

{Ct } = [<h ,y > ]{Ct 
-HT ti m n y 

(7.35) 


(7.36) 


(7.37) 


The formulations (7. 34) - (7. 37) have the advantage that the particular deflection 
basis can be factored from the solution process. In addition, the correspond- 

ing vector of downwashes may be computed by transforming the deflection basis to 
the downwash polynomials using (5.35) 

(:^ + ik){Y }= (^ + ik) [g'^ ]{a } = [g”^ ] (;^ + ik){a } 
dx • m dx mn n mn dx n 


and then differentiating using (5.40j) to obtain 


(:^+ik){y } = [G^ ] [d“ +ikS ]{a }, 

dx m mn mn mn n 


(7.38) 


or equivalently, by differentiating first using (5.402) 

(^+ik){y } = [D +ik6 ] {y } 
dx m mn mn n 

and then transforming to the downwash basis functions to obtain 


(— +ik){y } = [dY + ikS ] [G^ ]{a }• 
dx — — — — - 


mn mn n 


(7.39) 


Since all matrices in (7.38) and (7.39) are lower triangular, the resulting matrix 
product 

[W ] = [gY ] [oa + ik6 ] = [d’^ + ik(S ] [g’^ ] 


mn mn 


mn mn 


(7.40) 
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must be lower triangular also. Following (7.38), 


[W ] 
mn 


= 2 


Following (7.40), 


[W ] = 2 
mn 
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One o'ntains the same result both ways. 


[Wmn^ = 2 


ik 

0 

0 

0 

l+2ik 

ik 

0 

0 

3+2ik 

2+2ik 

ik 

0 

6+2 ik 

5+2ik 

3+2 ik 

ik 

10+2ik 

9+2 ik 

7+2ik 

4+2ik 


0 

0 

0 

0 

ik 


(7.41) 


The matrix [W ] represents the linear transformation of downwash from the deflec- 
mn 

tion basis {y } to the downwash basis {a }. 

m n 


(^+ik){y } = [W ]{o } 
dx m mn n 


(7.42) 


It is nonsingular if and only if the frequency is nonzero. 



Referring back to equation (7.21), we see that the aerodynamic work matrix 
corresponding to this basis is given by 


[a ] = -r . Y > ] ' 
mn 2 ■ m n Y 


(7.43) 


where ip satisfies Bland's integral equation for the basis function hjj^ = Y » 


(c— + ik) Y = T ip . 
dx m m 


(7.44) 


In the special case of steady flow in an infinite atmosphere, (5.13) states 

that the Fourier coefficients of pressure in the basis {y } must equal the Fourier 

n 

coefficients of downwash in the basis {a }; i.e., 

n 


<\p„,Y > 
^m' ‘n Y 


<(— + ik) 
dx 


a > 
n a 


W . 

mm 


This gives 

(k=0 & ^= 00 ) r^] = I (7,45) 

which reduces by inspection to equation (7.16) that previously was derived without 
benefit of the above matrix developments. 

In general the aerodynamic work matrix is nonsymmetric. Physically this is 
because the interaction between the air and the airplane structure is nonconservative. 
Mathematically it cannot be guaranteed that an arbitrary nonsymmetric matrix is con- 
gruent to a triangular matrix. In other words on purely mathematical grounds there 
need not exist a deflection basis in which the aerodynamic work matrix is triangu- 
lar. However it has been shown above in (7.16) and (7.45) that if the pressure poly- 
nomials are used as a deflection basis then the aerodynamic work matrix is in fact 
upper triangular for steady flow in an infinite atmosphere. At the present time, 
the meaning of this result is not fully clear, but it does cast the airfoil polynom- 
ials in a special role. 

A somewhat more abstract algebraic interpretation may be given to the transform- 
ation properties (7 . 29) - (7 . 32) of the unsteady airloads. In the case of linear 

aerodynamics, the pressure factor depends linearly upon the airfoil displacements. 

2 

Therefore, by the Riesz representation theorem, the vector space Ly of pressure 
factors is contained in the dual space L^* of airfoil displacements 
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In this regard the section coefficients are functions which map elements of the 
dual space linearly into the complex numbers 



* -*■ C 

* -y C 


and are therefore pure contravariant tensors of order 1 relative to [46,p.l2l. 

In the same regard, the aerodynamic work matrix is a function which maps pairs of 

2 2 

elements from the vector space of displacements L-^ and its dual space bilin- 
early into the complex numbers, i.e. , 

X i.2* ^ c. 


and is therefore a mi^ed tens or o f contravariant order 1 and covariant order 1 rela- 
tive to ^ [op cit] . 
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§8. Representation of airfoils with multiple controls 

In section 7 we discussed how unsteady airloads can be calculated for a general 
airfoil displacement in terms of the airloads corresponding to an arbitrary displace- 
ment basis and that certain advantages accrue if the airfoil pressure polynomials 
{Yin} are used as the canonical basis for displacements. This section describes how 
displacement functions {hm) possessing discontinuities of the type found with multi- 
ple leading and trailing edge controls can be specified by discrete input data, and 
how the matrix of Fourier coefficients 

[Hmn} = 

can be calculated for such functions. 


8.1 Input data specification of displacements 

Consider an airfoil displacement function of the type illustrated in Figure 1. 
Such a function may represent an upper, lower or mean profile. There may be multi- 
ple leading and trailing edge controls with sealed gaps. Discontinuities in the 
first derivative of the displacement function usually correspond to control deflec- 
tions. Discontinuities in higher derivatives correspond to changes in curvature due 
to design and/or aeroelastic effects such as changes in stiffness properties, etc. 
Values of displacement will be specified at distinct points 

X2 , . . . , 

which are called nodes . We shall assirme without loss of generality that nodes are 
labeled from left to right along the chord. 

Xi<...<Xn^ (8.1) 


The number of nodes may range from a minimum of one to some finite number , 

■^max 


1<N. 


<N 


^max' 


( 8 . 2 ) 


In addition, nodal values of different displacement functions 

hjj' • • ■ 

are to be given as input data. Some of the nodes may also be hinges at which dis- 
continuities in derivatives of displacement can occur. Hinge points must lie strictly 
between the leading and trailing edges. 
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-The total number of hinges may range from zero if there are none to at most 

Ng 

°max 


0<Ng£N 


iSmax 


( 8 . 3 ) 


( 8 . 4 ) 
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8.2 Subdivision of the chord into line elements 

This subsection describes the procedure for computing interpolation functions 
of airfoil displacements with hinge type discontinuities. The scheme employed 
subdivides the chord into line elements and is illustrated in Figure 4 for two 
hinges . 


k c=: 


e=2 


e=3 - ^ 



Figure 4. Subdivision into line elements 


The number of line elements along the chord equals one more than the number 
of hinges. 

Ng = 1 + Ng (8.5) 

If there is only one element, the element domain lies between the leading and 
trailing edges. 

SIq = {x : -l<x<l} if e = 1 & Ng = 1 (8.6) 

If there are more than one element, the first and last domains lie respectively 
between the leading edge and the first hinge, and between the last hinge and the 
trailing edge. 

fig = {x : -l<x<Xj^^} if e = 1 & Ng>l (8.7) 

= {x : Xn.< x<l} if e = Ng & Np>l 
c * N j c c 

u 


( 8 . 8 ) 
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All remaining domains are bounded by hinges. 

= (x : Xn <x<x_ > if Ke<N . (8.9) 

e ^-1 e 

For theoretical purposes, these domains are disjoint, i.e., 

f2ei ^B2 = <l> if ex e2, 

their closures intersect only at hinges , 

^ei ^e+1 ~ ^^ne^ ’ l^e^Ng— 1, 

and the interior of the union of their closures equals the entire chord between 
the leading and trailing edges. 

N _ 

(^=int( fig) = {x:-l<x<l}. 
e=l 

These properties are those required of finite element models [ 47 ,Ch. 6 ]. Thus, 
the interpolation procedures for the present method of solution in TWODI are con- 
sistent with the finite element method. The advantage of the element by element 
formulation above is that the interpolating functions for each element may be de- 
termined separately and by a uniform procedure. 

The numbers Nj.(e) of nodes contained in the closure of the various line 


elements are given by 

Njj(e) = if e = 1 & Ng = 1, (8.10) 

Nj^(e) = nx if e = 1 & Ng > 1, (6.11) 

Nx(e) = 1+n -n if Ke<Ng & N„ > 1 , (8.12) 

e e-1 ^ ^ 

Njj(e) = l+N^-n^_^ if e = Ng & Ng > 1. (8.13) 

The particular nodes x^(e) which are contained in fig are given by 

X (e) = x_ if l<n<N^(e) & e = 1, (8.14) 

n n ^ 

Xf^(e) = ^ if l<n<N^(e) & 2<e£Ng. (8.15) 


In this manner, the last node of one element is the same as the first node of the 
next element, as is indicated in Figure 4. 

The nodal values of deflections within each element are given in similar fashion. 

h (e) = h if l<m<N, & l<n<N (e) & e = 1, (8.16) 

mn mn h x 


h (e) 
mn 


h 


m, n-l+n 


e -1 


if Km<N, 
h 


& l<n<N (e) 

X 


2<e<N . 
e 


(8.17) 
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The interpolation functions which will be employed to represent displacements 
are airfoil polynomial splines 


Nx(e) 

h®(x) = I Y^(x), l£m<Nh, l£e£N^. (8.18) 

At each element nodal point x^(e)efi^, we require that the value of each inter- 
polation function equal the value of displacement prescribed by input data. Thus, 


h 

mn 


(e) 


N (e) 

X ^ 

y Yo(x (e))H „(e), Kn<N (e) , 
a n mJl X 


l<m<R , Ke<N . 

n ^ — e 


(8.19) 


Since the airfoil polynomials are unisolvent [5,p.42] each coefficient matrix 

[Yj(x^(e))l, l££,n£N^(e), l£e£N^ 

possesses an inverse, and since the nodal values, 

h (e) , l<m<N, , l<n<N (e) , l<e<N , 

are given by (8.16) and (8.17), the coefficients 

H „ (e) , l<m<N, , 1<£<N (e) , Ke<N 
m£ h X e 

are computable quantities. 
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8.3 Use of jump functions 

Consider the displacement interpolation functions for two adjacent elements: 


h (x) = y H (e) y.{x) if X E , 

T" ^ mZ i, “ 


Z=1 


e+i Nx(e+1) ^ 

h (x) = y H (e+l)Y (x) if X e fJ . 
m i e+1 


The expressions above are not yet in a form suitable for the solution process of 
TWODI because they are not uniform over the entire chord. This can be accomplished 
by utilizing jump functions and leads to infinite Fourier series whenever a discon- 
tinuity exists in a derivative. 

A representation for h which is uniformly valid over both line elements 

m 

X e int (SI u , ) 
e+1 


is given by 


N (e) 

X 


n=l 


H (e)v (x) + 
mn n 


N (e,e+l) N (e) H (e+l)-H (e) 

^ X mi ma (k) , 

Yd 


k=l 


y 

i=k+l 


kl 


) <x-x 
n n 

e 6 


where 


N^(e,e+1) = max{N^(e), N^(e+1)}-1. 


( 8 . 20 ) 


The function above equals h^(x) .and all its derivatives whenever x e S) and 

e+1 ® 

equals h^ (x) and all its derivatives whenever x e Slg+l- At the hinge x = ^m^ » 

the function is continuous in accordance with the sealed gap condition, and will 

e 

possess jumps in derivatives exactly equal to jumps in the derivatives from h^ to 


Fourier series expansions have been computed for the jump functions 

00 

k V k 


<x-a> = y <<x-a> ,7 > y (x) 


n=l 


n Y n 


( 8 . 21 ) 


as discussed above in section 5. In order to utilize these expansions, we extend 
the uniformization above to all chord elements by starting with the leading edge 
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element and adding jump contributions from all hinges dovmstream. It is easy to 
see that the result is 

Njj(l) 

h (x) = y H (1) Y (X) 
m mn n 

n=l 


“fi N (e,e+l) N„(e) H .(e+l)-H .(e) ,,, 

V ^ V mS- m^. (k) , V , .k 

+ A I Z yo )<x-x^ > . 


e=l k=l £.=k+l 


"•e ""e 


( 8 . 22 ) 


Equation (8.22) is uniformly valid for all points x e between the leading and 
trailing edges. 

(k) 

The derivatives Y„ (xn ) the hinges may be computed with the aid of equation 
X, -e 

(5.40) according to 


i-1 

£ni 



I 

ni=l 


(8.23i) 

£-1 

£-2 

“‘"1 °”l”2 ’'”2' 


I 

ni=l 

I 

n2=l 

(8.232) 

£-1 

I 

ni=l 

£-2 

I 

n2=l 

£-k 

- - • y Y 

n =1 ^"^1 "l‘^2 ^k-l^k rik, k>3. 

k 

( 8 . 233 ) 


Upon combining (8. 21) - (8. 23) , there results 


h (x) = y H Y (x) 

y-t bj mr*i 


(8.24) 


n=l 


where H is as introduced in (7.33) and is given by 
mn 


N, 


6 N^(e,e+1) N^(e) (etl) (e) 


H =H (1)+ I ^ I I 

e=l k=l Jl=k+1 


kl. 


(X )<<x 




(8.25) 


n Y 
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§9. Instructions for the use of TWODI 

This section describes the preparation of input data for use of the current 
TWODI program (TWODl-III) . it is assumed that the user knows the physical meaning 
of all terms used in this report and the procedures for accessing TWODI at his 
computer facility. Also, we reiterate that the solution is based upon the mathe- 
matical assumptions of inviscid subsonic linearized potential flow about a thin 
airfoil located midway between two parallel wind tunnel walls. 

9.1 Summary of capabilities 

TWODI will operate in either TIMESHARE (remote interactive terminals) or BATCH 
(noninteractive card jobs) . Input and output are fully compatible between TIMESHARE 
and BATCH. In addition TWODI is coded in ANSI FORTRAN to achieve maximum machine 
independence . 

TWODI will predict unsteady airloads consisting of any combination of the fol- 
lowing output quantities: 

(1) Fourier coefficients of pressure, 

(2) Values of pressure, 

(3) Section coefficients and center of pressure, and 

(4) Aerodynamic work matrix. 

The primary parameters determining the standard solution output are: 

(1) Mach number, 

(2) Reduced frequency, 

(3) Height to chord ratio, 

(4) Wall ventilation coefficient, and 

(5) Airfoil profiles (i.e., displacement mode shapes). 

The solution process in TWODI handles one case at a time, defined by a single com- 
bination of Mach number, reduced frequency, height to chord ratio and ventilation 
coefficient. Since the downwash corresponding to any given mode shape can be factor- 
ed out of the solution equations, multiple downwashes are handled simultaneously 
within a given case. 

Multiple cases may be handled under a single problem which consists of all 


combinations of various numbers of Mach number, reduced frequency, height to 
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chord ratio and ventilation coefficient. Within a given problem, selection of 
output quantities remains fixed. 

Provision is also made to enter several different problems. This may be 
done by altering a previously defined problem or by entering a completely new 
problem. All data, including the selection of output quantities, may be changed 
from one problem to another, with the exception that all problems in a given run 
must be done in TIMESHARE only, or in BATCH only. 

It is possible to solve the same problem using more than one method. This 
provision is preparatory to allowing TWODI to select automatically the optimum 
method of solution and thus to handle its own accuracy control at some time in 
the future. In addition certain special purpose and checkout calculations are 
available. 

All input data are automatically checked for correctness. If the run is a 
BATCH job, unacceptable data are selectively deleted with an explanatory comment 
and execution continues to the extent that it can. If the run is an interactive 
TIMESHARE job the user will be prompted to correct unacceptable data. 
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9.2 Interactive input format for TIMESHARE 

Queries and messages by TWODI are denoted with Q, responses by the user are 
denoted with A. Consecutive data entries must be separated by a comma or by one or 
more blanks. All data are assigned initial values by TWODI using subroutine INITLZ. 
These initial values are indicated within brackets [ ] below and are precisely de- 
fined in the glossary in section 9.4. After one or more problems have been entered, 
the most recently defined values are those of the most recently entered problem. In 
general, for default to the most recently defined values, type D followed by an im- 
mediate carriage return. 

Data module 0. Introduction 
Q: This is TWODI-III 

FOR DEFAULT TO INITIAL OR MOST RECENTLY DEFINED VALUES 
TYPE D FOLLOWED BY CARRIAGE RETURN IF IN TIMESHARE AND ENTER 
AN OTHERWISE BLANK CARD WITH A D IN COLUMN 1 IF IN BATCH. 

IF IN TIMESHARE TYPE HALT TO STOP. 

Data module I. Run parameters 
Q: ARE YOU IN TIMESHARE OR BATCH? 

A: TIMESHARE [TIMESHARE] 

Q; ENTER NUMBER OF LINES PER PAGE 

A: Type integer number of lines per page or D followed by carriage 

return for default. [66] 

Data module II. Output parameters 
Q: ENTER TITLE 

A: Type descriptive alphanumeric title of 1 to 72 characters or D followed 

by carriage return for default. [SAT^PLE PROBLEM] 

Q: ENTER DESIRED OUTPUT COMBINATION OF FOURIER, SECTION AND WORK 

A; Type desired combination or D followed by carriage return for 
default. [FOURIER, SECTION] 

Q; ENTER LIST OF PRESSURE POINTS 

A. Use standard list format (refer to glossary in section 9.4). Type 0 if 
none or D followed by carriage return for default. [10/-. 8,1] 
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Data module III. Flow parameters 
Q: ENTER LIST OF MACH NUMBERS 

A: Use Standard list format or type D followed by carriage return for 

default. [1,0] 

Q: ENTER LIST OF FREQUENCIES 

A: Use standard list format or type D followed by carriage return for 

default. [2,0,1] 

Q: ENTER LIST OF HEIGHT TO CHORD RATIOS 

A: Use standard list format or type D followed by carriage return for 

default. [1, INFINITY] 

Q: ENTER LIST OF MASS EFFECT VENTILATION COEFFICIENTS 

A: Use standard list format or type D followed by carriage return for 

default. [1, INFINITY] 

Q: ENTER LIST OF VISCOUS EFFECT VENTILATION COEFFICIENTS 

A: Use Standard list format or type D followed by carriage return for 

default. [1,0] 

Data module IV. Modal parameters 
Q; ENTER LIST OF NODES 

A: Use standard list format or type D followed by carriage return for 

default. [3/-l,l] 

Q: ENTER NUMBER OF MODE SHAPES 

A: Type integer number of mode shapes or D followed by carriage return for 

default. [3] 

A: ENTER MODE SHAPE 1 

A: Type values of mode shape 1 at modal collcoation points or D followed 

by carriage return for default. [l//n”, IZ/tT , IZ/tT] 

Repeat the last Q&A for the remaining mode shapes. IZ*^, 3/yir] , [IZ*^, -IZ*^, 5Z»^] 

Q; ENTER LIST OF NODE NUMBERS OF HINGES 

A: Use standard list format if downwash discontinuities are present. 

Type 0 if there are none or D followed by carriage return for default. [0] 

Data module V. Method parameters 
Q: ENTER NUMBER OF METHODS OF SOLUTION 

A: Type integer number of methods of solution or type D followed by carriage 

return for default. [1] 
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Q: ENTER SOLUTION PARAMETERS FOR METHOD 1 

A: Type the parameters for the first method of solution or type 

D followed by carriage return for default. Refer to the glossary for 
precise definitions. [3, 5, 0,0,0] 

Repeat the last Q&A for the remaining methods of solution. 

Data module IV. Data editing 

Q: DO YOU WANT THE INPUT DATA LISTED? 

A: Type YES or NO. 

Q; DO YOU WANT TO MAKE CHANGES? 

A: Type YES or NO. 

Note: An answer of YES will result in the following sequence: 

Q; DO YOU WANT TO LINE EDIT? 

A: Type YES or NO. 

Note: An answer of NO will result in data modules II-VI being repeated. 

For each data item which does not need to be corrected, simply 
type D followed by carriage return. An answer of YES will result 
in the following sequence ; 

Q: NOW OPEN FOR LINE EDITING. WHEN DONE TYPE END. 

A: Enter keyword for data item you wish to change aind you will 

be prompted for the relevant input. The keywords are as 
follows : 

TITLE 

OUTPUT 

PRESSURE 

MACH 

FREQUENCY 

HEIGHT 

MASS 

VISCOUS 

NODE 

MODE 

HINGE 

METHOD 

END 

The first three letters of any keyword are sufficient. This 
process may be repeated as often as desired. To terminate 
editing enter the keyword END. 



-58- 


Data module VII. Multiple problems 

Q: DO YOU WANT TO ENTER ANOTHER PROBLEM? 

A: Type YES or NO. 

Note: An answer of NO will result in the current list of problems being run. 

An answer of YES will result in the query 

Q: IF YOU WANT TO MODIFY AN OLD PROBLEM, ENTER ITS NUMBER 

OTHERWISE TYPE D FOLLOWED BY CARRIAGE RETURN 
A: Type N to alter problem N. Type D followed by carriage return to 

begin a completely new problem. 

Once the last Q is answered with NO the list of problems will be run. Upon 

completion, prompting will continue as follows. 

Q: DO YOU WANT TO ENTER ANOTHER PROBLEM? 

A: Type YES or NO 

Note: An answer of NO will stop the program. An answer of YES will result 

in the query 

Q: IF YOU WANT TO RETAIN OR MODIFY PROBLEMS FROM THE OLD PROBLEM LIST 

ENTER PROBLEM N OF THE OLD LIST AS PROBLEM M OF THE NEW LIST WHERE 
M IS LESS THAN OR EQUAL TO N. OTHERWISE PROBLEM N WILL BE DESTROYED. 
TO MODIFY AN OLD PROBLEM, ENTER ITS NUMBER 
OTHERWISE TYPE D FOLLOWED BY CARRIAGE RETURN 
A: Type N to alter or retain problem N as problem 1 of the new problem 

list. Type D followed by carriage return to begin a completely new 
problem list. Up to three problems may be defined. 

The APPENDIX presents a sample interactive input/output. 
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9 . 3 Card input format for BATCH 

Data are read in free format. Entries on a given card may be arbitrarily 

spaced so long as correct order is maintained. Consecutive entries on a card 

must be separated by a comma or one or more blanks. Otherwise blank cards with 

D in column 1 imply default to the most recently assigned values. 

Data module I. Run parameters 

Card 1 : BATCH 

Card 2: Number of lines per page or otherwise blank card with D in 

column 1 for default value of 66. 

Data module II. Output parameters 

Card 1: Descriptive alphanumeric title in columns 1-72. 

Card 2: Desired combination of the words FOURIER, SECTION and WORK 

separated by commas or blanks. Use an otherwise blank card 
with D in column 1 for default. 

Card(s) 3: Pressure points in standard list format (see glossary in 

section 9.4). Enter 0 if no pressure values are desired. Use 
an otherwise blank card with D in column 1 for default. 

Data module III. Flow parameters 

Card(s) 1; Mach numbers in standard list format. Use otherwise blank 
card with D in column 1 for default. 

Card(s) 2: Reduced frequencies in standard list format. Use otherwise 

blank card with D in column 1 for default. 

Card(s) 3; Height to chord ratios in standard list format. Use otherwise 
blank card with D in column 1 for default. 

Card(s) 4: Mass effect ventilation coefficients in standard list format. 

Use otherwise blank card with D in column 1 for default. 

Card(s) 5: Viscous effect ventilation coefficients in standard list fo 2 rmat. 

Use otherwise blank card with D in column 1 for default. 

Data module IV. Modal parameters 

Card(s) 1: Nodes in standard list format. Use otherwise blank card 

with D in column 1 for default. 

Card(s) 2: Number of modes, or otherwise blank card with D in column 1 for 


default. 
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Card (s) 3: Values of mode shape 1 at nodes, or otherwise blank Ccird with 

D in coltamn 1 for default. 

Repeat this card(s) for the remaining mode shapes. 

Card(s) 4: Node numbers of hinges if such discontinuities are present. 

Enter 0 if none are present, or use an otherwise blank card with 
D in column 1 for default. 

Data module V. Method parameters 

Card 1: Number of different methods of solution, or otherwise blank 

card with D in column 1 for default. 

Card 2: Solution parameters for first method of solution, or otherwise 

blank card with D in column 1 for default. See glossary for precise 
definitions of names and parameters. 

Repeat Card 2 for the remaining methods of solution. 

Data module VI. Data editing 

Data editing is not currently permitted with BATCH input. 

Data module vii. Multiple problems 

Card 1 : YES or NO beginning in coliimn 1 . If YES then repeat data modules 

II-VII. If NO processing will terminate. 
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9.4 Glossary of input terminology 

This subsection provides in alphabetical order definitions of input parameters 
needed to run TWODI. All data are automatically tested upon input and unacceptable 
values are deleted with an explanatory message. 

BATCH. Indicates noninteractive jobs with punched card input. 

FOURIER. Entering the word FOURIER causes the pressure Fourier coefficients to be 
printed upon output. Initial default condition is affirmative. 

LIST FORMAT. There are two standard formats for entering lists of numbers, one 
for arbitrarily spaced data and another for equally spaced real data. To enter N 
arbitrarily spaced real or integer data with values VI, ...,VN, type 

N, VI, . . - ,VN 

Such an arbitrarily spaced list may contain one (and at most one) infinite value 
by entering INFINITY. To enter N equally spaced real numbers from a to b, type 

N / A, B 

To enter an empty list, type 0. 

LIST OF FREQUENC#ES. From 1 to 10 finite real numbers representing reduced fre- 
quency based on semichord. May not be an empty list. Initial default condition 
is two frequencies with values 0 and 1. 

LIST OF height TO CHORD RATIOS. From 1 to 5 positive numbers representing height 
to chord ratio. May not be an empty list. To specify free air conditions, type 
INFINITY. Initial default condition is one value of INFINITY. 

LIST OF MACH NUMBERS. From 1 to 10 values of subsonic Mach number, each of which 
must be non-negative and less than one. May not be an empty list. Initial default 
condition is one Mach number with value 0. 

LIST OF NODES. From 1 to 20 distinct points at which mode shapes are to be collo- 
cated. May not be an empty list. Initial default condition is 3 nodes with values 
-1, 0 and 1. 

LIST OF NODE NUMBERS OF HINGES. From 0 to 3 integers locating possible flap hinges 
at the corresponding nodes. May be an empty list. Initial default condition is 
that there are no hinges. 

LIST OF PRESSURE POINTS. From 0 to 100 points along the chord where pressures are 
to be printed upon output. Each value must be greater than -1 (i.e., aft of lead- 
ing edge) and less than or equal to +l(i.e., not aft of trailing edge). Should 
not coincide with any hinge locations; if so, such pressure points will be deleted 
and the number of pressure points reduced accordingly. May be an empty list. Ini- 
tail default condition is 10 points equally spaced from -.8 to 1. 
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LIST OF MASS EFFECT VENTILATION COEFFICIENTS. From 1 to 5 values of mass effect 
(slotted wall) ventilation coefficients. Enter INFINITY for closed walls. For 
open jet enter 0 for both mass effect and viscous effect (see below) ventilation 
coefficients. Initial default condition is one value of INFINITY. 

LIST OF VISCOUS EFFECT VENTILATION COEFFICIENTS.. From 1 to 5 values of viscous 
effect (porous wall) ventilation coefficients. Enter INFINITY for closed walls. 

For open jet enter 0 for both mass effect (above) and viscous effect ventilation 
coefficients. Initial default condition is one value of 0. 

NUMBER OF LINES PER PAGE. Indicates number of lines per page of printed output. 

Not presently utilized. Initial default value is 66 in both BATCH and TIMESHARE. 

NUMBER OF METHODS OF SOLUTION. Integer from 1 to 10 which controls the number of 
different methods by which each case is solved. Presently used for convergence 
and extrapolation control. Initial default value is 1. 

NUMBER OF MODE SHARES. Integer from 1 to 5 which controls the number of different 
mode shapes for which airloads are calculated. Initial default value is 3. 

SECTION. Entering the word SECTION causes section coefficients C-j^ and Cj^ for lift 
and moment and chordwise center of pressure XQp to be printed upon output. Initial 
default condition is affirmative. 

SOLUTION PARAMETERS.* Three integers II, 12 and 13 followed by two real numbers 
HI and R2. Used to control the method of solution and to perform special purpose 
and check calculations. Will be phased out in the future and replaced by automatic 
accuracy control. The five-tuple (II, I2 , 13, Rl, R2) produces the following (a 
dot • indicates parameters not presently used, so enter 0) : 

(1,2,1, •,*) Checks the Legendre-Gaussian quadrature tables. 

(1,2,2, •,■) Checks the logarithmic-Gaussian quadrature tables. 

(1,2,3, •,•) Checks the inverse square root Gaussian quadrature tables. 

(1,2,4, •,•) Checks the Laguerre-Gaussian quadrature tables. 

(1,3, •,•,•) Checks the eigensolution of ( 3 .II 4 ). 

(1,4, •,*,*) Checks the infinite series summation (3.11j). 

(1,5, •,•,•) Checks the calculation of K^; in (3.10). 

(1,6, •,*,•) Checks the calculation of Fj in (12.26). 

(1,7, •,*,•) Checks the calculation of F 2 in (12.33). 

*Note added in proof. Since the writing of this report, TWODI has been modified 
for automatic extrapolation as described in section 11 and the Possio free air 
kernel is computed as described in section 12. This Glossary, the Addendum 
immediately following the Conclusions in section 15 and the Appendix reflect these 
changes . 
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Checks the calculation of airfoil polynomials in (5.5). 

Checks the calculation of the collocation matrix (6.8). 

Special purpose. Calculates the Kiissner-Schwarz solution for 
flaps with hinge at Rl and reduced frequency R2. -1SR1<1. 

Special purpose. Quadrature evaluation of Possio kernel. 
Special purpose. Quadrature of x log x using Legendre-Gaussian 
quadrature . 

Standard solution. Richardson extrapolation of order 13 
based on Bland ' s collocation method with a period of 12 
basis elements using internal error tolerances of Rl. If 
R1=0, Rl is set to 10~^. Final solution error is approxi- 
mately lORl. Reduces to Bland's collocation method with I2 
basis elements when 13=0. 1<I1<24. This is the only combi- 

nation of solution parameters needed to compute airloads. 

Use of solution parameters with II 3 requires expert interpretation and is not 
recommended for production usage. Initial default value of solution parameters 
is (3, 5, 0,0,0). 

TIMESHARE. Indicates remote terminal jobs in which TWODI and the user communicate 
interactively. Initial default value is affirmative. 

WORK. Entering the word WORK causes the aerodynamic work matrix [Amn] to be 
printed upon output. Not presently implemented. Initial default condition is 
negative. 


(1,9, •,•,•) 
( 1 , 10 , •,•,•) 
(2,1, •,R1,R2) 

( 2 , 2 , •,•,•) 

( 2 , 4 , •,•,•) 

(3,I2,I3,R1, •) 
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§10. Convergence character istics, of T WODI 

In our previous report [5] an extensive examination of the convergence 
of Bland's collocation method was made; both analytically and computationally. 

Our principal theoretical result was that collocation and Galerkin’s method were 
numerically equivalent, so that convergence was analyzed in terms of Galerkin's 
method. On the assumption that the quadrature errors in the evaluation of the 
Galerkin matrix could be neglected, it was concluded that Bland's method was con- 
vergent. This theoretical result was validated by exhibiting the rapid numerical 
convergence of TWODI-I for a variety of steady and unsteady flow problems with 
smooth downwashes. 


In the present section these results are extended in several directions. 

First our convergence result is strengthened by showing that one need not make 
the assumption of negligibility of quadrature errors in the Galerkin matrix and 
we present a direct proof of the convergence of collocation. This enables us to 
obtain improved error estimates along with strenthening the theoretical basis of 
collocation methods in aerodynamic calculations. Second, the numerical conver- 
gence of TWODI is demonstrated for problems with flaps. An important result in 
this case is the verification of the predicted slow rate of convergence, in contrast 
to the rapid rate achieved for smooth downwashes. This observation required us to 
make a detailed study of a number of convergence accelerating techniques; a topic 
which is taken up in detail in the following section. 


10.1 l 2 convergence of collocation. 

Since we wish our analysis to cover other equations than those with Bland's 
kernel we begin by making several assumptions concerning the kernel K (x) =K (x) -1/x. 
These are: 

(A-1) K(x- 5) is Lebesgue-Stieltjes square integrable with respect to the product 




Ogi / 1+X 


/ 1 +^ 


|k(x- 5) |2 dCdx<=». 


(A-2) K(x-^) is mean square continuous,- i.e.. 



K(x-C+h)-ic(x-C) I 2 dS 


0 . 


( 10 . 1 ) 


lim 

h^ 


( 10 . 2 ) 
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For the next assumption we need the following definition. 

Definition . Let feL^ and let Qjj(f) denote the sxmi given by the quadrature 
rule (5.30^). We say that f is quadrature convergent for tQjg} if 
defined, N^l, and if 


lim QN(f) 
N-x® 


i: 



f (x)dx. 


(A-3) For each fixed |k^Cx)| = |K(x-g) | is quadrature convergent. 

Since for Bland's kernel K(x) = (x) + K (x) where K (x) is continuous 

L c c 

it follows by standard theorems [44] , [48] that satisfies (A-l)-(A-3). For 

Kjj (x) 1 similar properties may be shown by direct, though lengthy computation 
[ 49 ]. The proof, particularly of (A-3), relies on the fact that Qjj(f) may be 
viewed as a Riemann sum [49], [50] and that appropriately chosen Riemann sums con- 
verge to the integral of (log|x|)^. 

Our next step is to recast the collocation method in a suitable abstract 


form. 

Definition . Let feL^ and let the zeros of Then 

N 

Ljj(f) = I Jlj^(x)f(x ), (10.3) 

k=l 

is the unique polynomial which interpolates to f (x) at Here 


are the fundamental polynomials of Lagrange interpolation. 

Using (10.3) the sequence of operators {Kj^} is defined by 




(10.4) 


Let 


rp (F) 


N 

= I 

k=l 




be the collocation approximation to ij) using theN basis elements 
a little algebra shows that satisfies the equation 


Then 


H'l'.. 


+ K ll) = L (w) — W„. 
N^N N N 


(10.5) 


^For technical reasons we shall consider log|x-Cl - 0 if x = ?• 
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(Note that (10.5) is analogous to the equation satisfied by iIj , where ip 

N N 

is the corresponding Galerkin approximation to ip.) The convergence of ijj to 

N 

ip is based on the following lemmas. 

Lemma 10.1 Assume that {Kjj} are bounded and that 

lira II K-K 11 = 0. (10.6) 

N 

N-^ '■ 

Then there exists an N such that for all N>N , (H+K ) exists, is bounded and 

o o N 

the norms = 1|(H+K^)”l|| are uniformly bounded. From this it follows that 

for all N>N , \p exists and 
o N 

II a- 

Proof . The existence and boundedness of (H+K^) ^ follow along standard lines 
[51]. To obtain (10.7) observe that 


^ =.|;-(H+K^) 



(H+IS^)"'[(H+KJ^)'^-WJ^] 


-1 ■ 

= (H+K„) [HI|J+L^, (KlJ>) -L^, (w) ] 

N N N 




[Hlp-L^(Hi;)) ] . 


( 10 . 8 ) 


Taking norms on both sides of (10.8) gives (10.7) □ 

Lemma 10.2. Let f^L^ and assume that |fp is quadrature convergent for 
Then 


lim 

N->-» 




a 


0 . 


(10.9) 


Proof . In [ 44 ] it was shown that the lemma holds when f is assumed to Riemann 
Stieltjes square integrable, in particular when f is continuous. However, a 
careful examination of the proof presented there shows it is sufficient to re- 
quire that I f I ^ be quadrature convergentD 

Using Lemmas 10.1 and 10.2 it follows that provided that [|k-K^|[ 0, 

K satisfies (A-1) - (A-3) and |wp is quadrature convergent for 
To see this, observe that it follows from the above discussion that \}i}p\^ is 
quadrature convergent since Hij; = Kij; - w and | Ki[i | ^ is quadrature convergent (in 
fact continuous from (A-2)) and lw|^ is quadrature convergent by assumption. 

Thus by Lemma 10.2 

2 

lira II Hi|)-L^(Hii;) II = 0^ (10.10) 

N-xx) 
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and so (10.7) shows that in the norm of L^. 

N Y 

our main convergence theorem. 

Theorem 10.1 . Let K(x-C) satisfy (A-1) - (A-3) . 
ture convergent , then 


Using this we arrive at 
Assume that w is guadra- 


lim II 'l'-'l'„|| = 0. 

Proof. From the preceeding discussion it suffices to show that are 

bounded and that 


lim II K-K^ll = 0. 


N->«o 


The boundedness of is established first. 

From (10.3) it is seen that „ 

N 

(K^jp) (x) = J ilj^(x) (Kifi) (X]^) 
k=l 

From Eq. 14.2.4 in [ 44 ] we get that 


|L^(f) fd6(x) = K 


ir 

where f “ T — ^ are the weights in the quadrature rule Q . Thus 


de(x) = I Xj^|ki|;(Xj^) 1^. 

-1 


But 


1kiJ;(x^) 1 - 


lK(x^-C)i|;(?)da(5) 


< ( 


|K(x^-C) I dct(4) ) ( 


-1 


I < 1^(0 I da(C)) , 


-1 


so that 


2 N 

■■ ( I 

k=l 




“In what follows we use da(£) = d6(x) - y y— 


dx. 
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where 


This gives 


,1 




-a 


|k(Xj^-?) I da(C). 




since ^ and so is bounded. 


To prove the uniform convergence of K to K observe that 

N 






where 


Thus 


N 

= K(x-?) - I £ (x)K(x -5), 
N , a k k 

k=l 


'-\ll i 


,1 ,1 


■1 -1 


|h^(x, 5)| d3(x)da(C). 


Let y^(x) = K(x-?). Then 


= y^(x) - L^(^J^)(x). 


From Lemma 10. 2 it follows that 

,1 


lim 

N->-co 




I H (x, C)| ^ dB (x) = lim 
N J ^ 


ly^-L^(y^) I d3(x) = 0 


for each fixed C- 

From the proof of Lemma 10.2 found in [44] we get 

,2 

. ^ V u 'I ^ Al 


2 f .. 2 

lim sup([ |y;--L (y )| d3(x))< 4 |K(x~^)j d3(x). 

J 5 N C - J 

1 -1 


Since 


1 .1 


|ic(x-?) I d3(x)da(C) 


-1 -1 


( 10 . 11 ) 


< Oo 
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it follows from {10.11) that there exists a fimction p(?) such that 



(x,5) 1 dB(x) £p (5) 


and 



so that by the dominated convergence theorem 


lim 





2 

d6(x)da(C) 


[lim 


(I |h^(x,5) I dB(x) )]da(e) 
-1 


0 . 


Thus 


lim 


K-K . 


N " — 


.1 .1 


lim 

N-hjo 


[H^(x,?) 


- 1-1 


2 

d6(x)da(C) 


0 


and the theorem is p roved D 

We now make several observations concerning the convergence theorem. First 
we note that it has been shown that converges to \p in mean square and not 
pointwise, although numerically pointwise convergence is indicated [ 5 ] - How- 
ever it is easily shown that the generalized Fourier coefficients converge 

to the true Fourier coefficients of i|). In addition integrated aerodynamic forces 
such as lift and pitching moment also converge to their true values. This was 
established in [ 5 ] and follows easily using the fact that such quantities are 

represented in terms of inner products of the pressure factor and an appro- 
2 

priate function in L . 
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Second we observe that the rate of convergence is proportional to 
]] lljjj • If fh® downwash is smooth, then in general Hi/; will be smooth 

and we expect L^(Hi{') to converge rapidly to Htp. This was observed in [5]. 

However if w is not smooth then Hi|; will be poorly behaved and thus slow con- 
vergence of to Hi/' is anticipated. In fact for the particular case of 

steady flow in free air, k=0 so that 

- II II “• (10.12) 

Although we have not obtained exact estimates of ||w-L^(w) ||^ we expect that 
asymptotically it should behave no worse than 0(1//N). The slowest conver- 
gence rate should occur for a downwash w(x) corresponding to a simple leading 
or trailing edge flap. Numerical results presented in Tables 3-4 indicate a 
somewhat better rate of 0(1/N). For the same rate is expected since 

II 'i^-^Nlly- ^nII II a- Vll II a ""ll !!« ^ 

and by (A-1) Kij; is continuous. Thus the dominant error term is t^|| w-L^ (w)| | . 

Again the results exhibited in Tables 5-6 indicate an 0(1/N) rate of convergence. 
In general for 24 basis elements, l%-2% error might be anticipated using 
TWODI for flaps. This would be an unacceptably large error for high precision 
engineering work. Substantial increase in accuracy thus appears to require one 
of two strategies; an increase in the number of basis elements and/or the 
utilization of alternate solution methods. 

Since present engineering technology allows pressure measurements to be 

-3 

made within an error of 0 (10 ) this is the maximum error that we would like to 
have. From Table 3 we see that an error of 0(10 for a midchord flap would re- 
quire several hundred basis elements. Such an increase in the number of basis 
elements is out of the question. The second option is pursued in the following 
section where it is shown that an error of 0(10“'^) can be achieved for a mid- 

_o 

chord flap using 16 basis elements, and an error of 0(10 ) is obtained for a 

three-quarter chord flap using 12 basis elements. 
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Table 3. Principal error term in norm for a inidchord flap (M=0 ,k=0 . n=°°) 


NP 

ll’J'jjpll Y 

% error 

1% 

error [ xNP 

1 

1.77245 

-10.545 


10.5 

2 

1.50774 

+5 . 964 


11.2 

3 

1.67441 

-4.430 


13,3 

4 

1.54963 

+3.352 


13.4 

5 

1.64833 

-2.843 


14,2 

6 

1.56602 

+2.330 


14.0 

7 

1.63626 

-2.051 


14.4 

8 

1.57475 

+1.785 


14.3 

9 

1.62929 

-1.617 


14.5 

10 

1.58017 

+1.447 


14.5 

11 

1.62477 

-1.335 


14.7 

12 

1.58386 

+1.217 


14.6 

13 

1.62159 

-1.136 


14.8 

14 

1.58654 

+1.049 


14.7 

15 

1.61923 

-0.989 


14.8 

16 

1.58857 

+0.923 


14.8 

17 

1.61471 

-0.876 


14.9 

18 

1.59817 

+0.823 


14.8 

19 

1. 61597 

-0.786 


14.9 

20 

1.59145 

+0.743 


14.9 


II II 7 exact 

/-IT 

= / 2 ^ 1 

= 1, 

.60337 


Table 4. Princioal 

error term 

in lift for 

a midchord flan (M=0,k=0.n=°°) 

NP 

Cl 

% error 

j % error | xNP 

1 

6. 28319 

+22. 203 

22.2 

2 

4.54656 

-11.573 

23.1 

3 

5.60728 

+9.057 

27.2 

4 

4.80272 

-6.591 

26.4 

5 

5.43401 

+5.687 

28.4 

6 

4. 90481 

-4.605 

27.6 

7 

5. 35469 

+4. 145 

29.2 

8 

4.95964 

-3.539 

28. 3 

9 

5. 30922 

+3.260 

29.3 

10 

4.99386 

-2.873 

28.7 

11 

5.27974 

+2.687 

29.6 

12 

5. 01725 

-2.418 

29.0 

13 

5. 25908 

+2. 285 

29.1 

14 

5.03424 

-2,088 

29.8 

15 

5.24379 

+1,988 

29.4 

16 

5.04715 

-1.837 

29.4 

17 

5.23202 

+1-759 

29.9 

18 

5.05729 

-1.640 

29.5 

19 

5. 22269 

+1. 577 

30.0 

20 

5.06546 

-1.481 

29.6 


C 


^exact 


2+if 
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Table 5. 

Cl for an 
6 

oscillating 

flap hinged at 

the 50% chord 

(M=0, k=.l,n=«>) 

source 

real 

imag 

magn 

phase 

% error (magn) 

NP=1 

5.24189 

-0.50914 

5.26656 

-5.55 

+20.9665 

NP=2 

3.83875 

-0.39049 

3.85856 

-5.81 

-11.3735 

NP=3 

4.71938 

-0.46530 

4.74226 

-5.63 

+ 8.9240 

NP=4 

4.05062 

-0.40762 

4.07108 

-5.75 

- 6.4922 

NP=5 

4.57548 

-0.45292 

4.59785 

-5.65 

+ 5.6071 

NP=6 

4.13538 

-0.41479 

4.15613 

-5.73 

- 4.5387 

NP=7 

4.50954 

-0.44719 

4.53166 

-5.66 

+ 4.0868 

NP=8 

4.18095 

-0.41868 

4.20186 

-5.72 

- 3.4883 

NP=9 

4.47172 

-0.44390 

4.49370 

-5.67 

+ 3.2149 

NP=10 

4.20940 

-0.42113 

4.23041 

-5.71 

- 2.8326 

NP=11 

4.44720 

-0.44177 

4.46909 

-5.67 

+ 2.6496 

NP=12 

4.22884 

-0.42281 

4.24993 

-5.71 

- 2.3842 

NP=13 

4.43001 

-0.44027 

4.45184 

-5.68 

+ 2.2534 

NP=14 

4. 24298 

-0.42403 

4.26411 

-5. 71 

- 2.0585 

NP=15 

4.41730 

-0.43916 

4.43907 

-5.68 

+ 1.9601 

NP=16 

4.25371 

-0.42495 

4. 27489 

-5.71 

- 1.8109 

NP=17 

4.40751 

-0.43831 

4.42925 

-5.68 

+ 1.7345 

NP=18 

4.26214 

-0.42568 

4.28335 

-5.70 

- 1.6166 

NP=19 

4.39974 

-0.43764 

4.42145 

-5.68 

+ 1.5554 

NP=20 

4. 26894 

-0.42627 

4.29017 

-5. 70 

- 1.4600 


Table 6. 

Cr for an 

oscillating 

flap hinged 

at the 75% chord 

(M=0 ,k=. 1 , ri=oo) 

source 

real 

imag 

magn 

phase 

% error (magn) 

NP=1 

5.20343 

-0.76932 

5.25999 

-8.41 

+62.2622 

NP=2 

3.80393 

-0. 57961 

3. 84783 

-8.66 

+18.6994 

NP=3 

2.86126 

-0.44037 

2.89495 

-8.75 

-10.6954 

NP=4 

4.01530 

-0.60740 

4.06098 

-8.60 

+25.2747 

NP=5 

3.45638 

-0,52720 

3.49635 

-8.67 

+ 7.8568 

NP=6 

3.01247 

-0.46186 

3.04767 

-8.72 

- 5.9843 

NP=7 

3.70555 

-0.56291 

3.74806 

-8.64 

+15.6216 

NP=8 

3.36402 

-0. 51342 

3.40298 

-8.68 

+ 4.9765 

NP=9 

3.07119 

-0.47035 

3.10700 

-8.71 

- 4.1540 

NP=10 

3.56680 

-0.54280 

3.60786 

-8.65 

+11.2967 

NP=11 

3.32122 

-0. 50707 

3.35971 

-8.68 

+ 3.6417 

NP=12 

3.10241 

-0.47489 

3.13854 

-8.70 

- 3.1811 

NP=13 

3.48818 

-0.53135 

3.52842 

-8.66 

+ 8.8461 

NP=14 

3.29653 

-0. 50341 

3.33475 

-8.68 

+ 2.8717 

NP=15 

3.12177 

-0.47771 

3,15811 

-8.70 

- 2.5774 

NP=16 

3.43758 

-0.52397 

3.47729 

-8.67 

+ 7.2688 

NP=17 

3.28046 

-0.50103 

3. 31850 

-8.68 

+ 2. 3704 

NP=18 

3.13496 

-0.47963 

3.17144 

-8.70 

- 2.1662 

NP=19 

3.40229 

-0.51882 

3.44162 

-8.67 

+ 6.1684 

NP=20 

3.26916 

-0.49936 

3. 30708 

-8.68 

+ 2.0181 
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10,2 Stability and integration error. 

The error estimate given by (10.7) is a theoretical one based on the assump- 
tion that all arithmetic is performed exactly. In general there will be other 
sources of error, particularly roundoff and errors in the evaluation of the col- 
location matrix. To see the effect of these assvmie that ijijg is the computed value 
of Then letting seen that 

(10.13) 

Taking norms on both sides of (10.13) shows that 


I 'J'-'/'N [1^1 II ’^'■\ll^ +11 '^’^'^jlly 


(10.14) 


Since || is estimated in (10.7) we concentrate on evaluating H'^’f^fjlly ’I'fj 


N "y 

is given by 


N 

= Jay 

n' n 

n=l 


(10.15) 


where {a } , solve the collocation equations (6.7), Let {a } , be the numeri- 

n m=l n n=l 

N 

cally computed values of {a } .. Then 

n n=l 


% = Vn- 


n=l 


This gives 


N ^ N • 

6 i(( = ^ (a -a )y = J Sa y . 

n n n n n n 

n=l n=l 


Using the orthonormality of 9®^ that 


N 2 i 

l%l!y= (I l«aj )2 


n=l 


(10.16) 


Letting 6 a = { 6 a^}|^_^, |16 i(|jj||y= || 5 aj |2 , where || 6 a ||2 is the usual Euclidean length 
of a complex N- vector. 


If C denotes the collocation matrix, then 
N 

Cn£= w 


(10.17) 
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and 


(C +AC )a = w, 
N N — — 


(10.18) 


where a = {a a = {a w = {w(x„)}^ and AC is the error in the evalu- 

— n n=l — n n=l — K K=1 N 

ation of C . Standard error estimates show that 
N 


l«5a]|^ <a|] AC^I 


-li 


(10.19) 


where o is independent of N [51] . Consequently the progagation of error is de- 
termined by turn to the problem of estimating this quantity. 

Rather than analyze ||c ^|j directly we introduce a suitable vector norm on 
C (C = {(z, , , . . . , z ) z. complex, i = 1,2,...,N}) and use the induced matrix 
norm instead. 

Definition. Let {A, , be the weights of the quadrature rule Q„. Let ze(E^ 

k k=l N 

and define 


izii® = <! 

k=i 


( 10 . 20 ) 


Since l]t>0 it is straightforward to verify that|| z|| ^ is a norm on (C^. Let 
T;C^-Kn'^ be an NxN complex matrix, and let |1t|| ® be the matrix norm of T induced 
by II II 2 ’ Using these definitions we arrive at the following theorem. 


Theorem 10.2. Let C_ be the collocation matrix and let t- siiP II (H+K ) 
N N - " N 

a=[a (x )], where {x, }, _ are the collocation points. Then 

n K k k=l 


-li 


Let 


|c '^11 2 < X 

' N " — 


-i|| Q 


Proof. Let a = (a and w = {w(x ) • Then 

— n n=l — K k=l 




Now 


so that 




n=l 


H(jj (x) = y a a (x) 

TJ “ nr) 


n=l 


and 


N 

(Ht]’ ) (x ) = y a a (x ) . 
N K n n K 

n=l 


( 10 . 21 ) 
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From (10.5) it follows that 


♦n - 


(10.22) 


so that 


II ’^Nll y ^ II (h+K^)"^|| ||v^!I^<t|| WjjII a = t( I I ^1 

' n=l 


I 2,1/2. 


Since H is unitary [ 5 ] ,|h|j„ | = , giving 

' N " Y ' N a 


Hi|) <T w 

I II w II 2 


Letting H = , (10.21) gives 

N Jx JC~" J- 


Taking norms we get 


But II h|^ = ||Hii) lla giving 


a = o H. 


I|a||« i ||c.-MI®IIh||§. 


y®lT||a-'|i°||w||fi 


(10.24) 


From (10.24) it follows that 


|c->|l®iT||a-l|l®D 


Using the theorem it is easily shown that 


Thus the propagation of numerical error in the collocation matrix depends essen- 
tially on IjcT^II At present we have no theory to predict the growth of 
llcT^II®. However numerical experimentation has shown that ||a~^||^ grows slowly 
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with N so that for 1<N<^20 the collocation matrix is well-conditioned. From 
this we conclude that reduction in numerical error for fixed N requires 

careful evaluation of the integrals in C . 

N 

In TWODT-I difficulties arose in unsteady problems where we found it 
difficult to obtain more than 3-4 decimal accuracy. Since the above analysis 
indicates that inverting the collocation matrix is numerically stable, we 
conjecture that the source of the difficulty is the improper integration of 
the log terms in the kernel. In section 12 we show that this is in fact the 
case. A thousand fold increase in accuracy is obtained by efficient integra- 
tion of the singular terms. 
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§11. Convergence Acceleration 


This section presents a theoretical error analysis of various methods which 
might be used for improving the rate of computational convergence when flaps are 
present. 

As we have shown in section 10 Bland's collocation method converges rapidly 
when the downwash is smooth. However, if the airfoil contains flaps, the results 
presented in Table 3-6 indicate that for physically important w's a convergence 
rate no better than o(l/N) can be expected. This is disappointing, since an en- 
gineering accuracy of 0(10“^) would appear to require something on the order of 
300 basis elements. As TWODI can presently accommodate a maximum of 24 basis 
elements it might appear that the flap problem is essentially intractable with 
today's computing capability. 

This circumstance naturally leads one to look for alternate methods of 
solving the genralized airfoil equation when flaps are present: the goal being 

to obtain a convergence rate sufficiently great so that the stated degree of 
accuracy can be achieved within the limitations of our existing code. 

A great many possibilities present themselves. Among the choices are: 

(a) Modification of the collocation method by either changing the basis elements, 
the collocation points, or both. 

(b) changing the basic method of solution to something like Galerkin's method 
[5] or least squares [8]. 

(c) Singularity subtraction. This is commonly referred to as Landahl's method 
[52] in three dimensional problems and has been developed by Rowe et al [53] , 
[54] . 

(d) Iterative improvement. 

(e) The use of reverse flow theorems [55] , [56] . 

(f) The use of extrapolation methods. 

Since Bland's collocation method is understood theoretically, and is efficient 
for smooth downwashes, it is presently felt that improvements should be sought 
utilizing as much of the output of TWODI as possible. For this reason we have 
ruled out category (a). In addition, our examination of results achieved by 
Milne in [56] indicates that the 0(1/N) rate of convergence of flaps may be in- 
herent in any collocation method using a continuous representation of regard- 
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less of the choice of colloation points or basis elements. However, some recent 
work by Nissim and Lottati [57] using a piecewise continuous representation of 
for Possio's equation seems to indicate that considerable improvement may be 
achieved over standard methods. Since we have had little time to evaluate this 
work, such representations are not considered in this report. 

Because of the three way equivalence of Galerkin's method, collocation and 
least squares established in [5] , use of methods in category (b) is also ruled 
out. We are thus left with techniques in (c)-(f). 

Although such methods have been discussed in the literature for many years, 
we have been unable to find any discussion of their theoretical properties and 
little evidence of the controlled numerical experimentation necessary to dis- 
tinguish among competing techniques. For this reason we have included a fairly 
detailed examination of these procedures which we hope will clarify our reasons 
for not finding them effective at the present time. Since the analysis that 
follows in sections 11.1 to 11.3 is fairly involved, we observe that the reader 
may, with little loss in continuity, skip these and proceed to section 11.4. The 
principal results reported in sections 11.1-11.3 indicate that the techniques 
(c)-(e) yield only a modest increase in the rate of convergence of Bland's colloca- 
tion method at the expense of substantially increased arithmetic. These facts 
ultimately lead us to examine the use of Richardson extrapolation as a means of 
accelerating convergence . 

While we have not yet completely automated this technique, it presently can 
be used in conjunction with some preliminary data analysis to obtain the stated 
goal of errors of order 10~^ or less for a wide variety of flow problems. The 
principal drawback to total automation is the oscillatory nature of the conver- 
gence when flaps are present. In addition, the period of oscillation appears to 
depend on the location of the hinge point, and cannot, at present, be predicted 
theoretically. As the analysis that follows shows, there is reason to believe 
that a combination of either singularity subtraction and/or iterative improvement 
in conjunction with extrapolation may yield the most efficient algorithms. 
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11.1 Landahl's method 

Because of its widespread use in three dimensional calculations [53] , [54] we 
begin with an examination of Landahl’s method for the solution of (6.1). The aim 
is to find an equivalent integral equation with a smoother downwash than w. In 
principal there are a variety of implementations of the basic idea, and we discuss 
two of these. From a mathematical point of view the method may be seen as a gen- 
eralization of the well known Kantorovich regularization method used for solving 
integral equations of the second kind [51] and our analysis follows closely that of 
Atkinson presented in [51] . 

11.1.1 Landahl I. 

We consider solving the integral equation 

(H+K) = w . (11.1) 

Let bs an approximation to ip. Using ip^ we define the residual downwash w^ by 

w = w-(H+K)t[) , (11.2) 

R o 

and the residual pressure by 

(H+K)i|< = w„. (11.3) 

K K. 


Theorem 11.1. tp = ijj +tp . 

Proof. From (11.2) and (11.3) it is seen that 


(H+K) (ill +iii ) = (H+K) ill +(H+K)ili = w-w +w = w. 
o R o R R R 

Thus \p +4> solves (11.1). Using the fact that (11.1) has a unique solution it fol- 
o R 

lows that ill = ill + ill n 
o R 

From Theorem (11.1) it is seen that i|i can be determined by solving (11.3) instead 

of (11.1). Basically we are considering a generalization of iterative refinement used 

for solving linear algebraic equations [51] . For the method to be effective ip^ must 

somehow be determined so that the residual downwash is smoother than w; thus solving 

for i|i would give a more rapid rate of convergence than obtaining [p directly. 

R 

Viewing the technique as iterative refinement we might proceed as follows. Let 

i|i be obtained using Bland's collocation with the maximal number of basis elements N. 
o 

(presently N=20) . That i , iji solves 


(H+K ) i[i 
N o 


w 


N‘ 


(11.4) 
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In this case 


w = (H+K) \b ’-w 
R o 


and letting if) be the approximation to if) using N basis elements gives 
R R 


(11.5) 


We take our approximation to if) as 


if) = if) +if) . 


( 11 . 6 ) 


Hence, (11.6) gives 




, -1 




= (H+Kj (L„(w)+L„((H+K)iii )-w)] 

N N N O 

= (H+K )~^[L (w)]. 

N N 

Thus if) =if) and no improvement can be expected. From the above result it is seen 
N o 

that if)^ must be chosen in some other way than by (11.4). 

To do this let H+K be decomposed as 


H+K = H+K^+K^ 

where it is assumed that (H+K^) ^ exists. Let 


(11.7) 


if)^ = (H+K^) ^w 


( 11 . 8 ) 


and assume that if)^ can be determined accurately (essentially analytically) . Then 
in this case 

w = w- (H+K) if) = w- (H+K, ) if) -K„if) = -K„iJj . 

R o lo.iO 2o 

Thus the residual pressure if) satisfies the eiguation 

R 


(H+K) if) = -K2if)g. 


(11.9) 


If it is assumed that K^ satisfies (A-2) , then w = -K_iJ) is continuous. If it 

2 R 2 o 

is possible to effect the kernel splitting K=K^+K 2 so that K^ is highly differen- 
tiable and (H+K^) has a known inverse, then (11.9) should present a more tractable 
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problem to solve than (11.1). To obtain a numerical approximation to if) (11.9) is 


solved by collocation giving a function ij) satisfying 

R 

■ -\<W - ”r 


( 11 . 10 ) 


Let 


ll; = iJ> +\b . 


(11.11) 


is then taken as our approximation to ijj. The following theorem justifies this 
procedure. 

Theorem 11.2. Let ip be defined by (11.11). Assume that K. satisfies (A-l)-(A-3) 

^ - 2 
so that II 0- Then converges to ip in L^ and 


11^111 (H+K^) ^||{|| 

i satisfies 
N 

- ”*<LnV L*0>- 


Proof. We first show that p satisfies 

N 


To see this observe that 


= HlP + Hil^ 

N R o 






- -SNV^lt.*E - 


= “K., b -K ill - K ll) +K, ill -K, p +w 
2N^N IN^R IN^o IN^o l^o 


( 11 . 12 ) 


(11.13) 


-K Ip + (K, -K, ) p +W. 
N^N IN 1 ^O 


Thus (11.13) follows. To obtain (11.12) write 


and 


i)j = (H+K) w 


iji„ = (H+K ) ^w+(H+K„)“^[(K^ -K)i{i ] 
N N N IN O 


(11.14) 
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Siabtracting gives 

= (H+K)“^w-(H+Kj^)"^w-(H+Kj^r^[(K^^-K^)l|;^] 

= (H+K^) (Kjj-K) (H+K) ~K- 

= (H+K^)“^{ (K^-K) (H+K)"^w-(K^^-K^) (H+K^)"^w} 

= (H+K )“^{(K -K) (H+K)"^-(K -K_ ) (H+K_)“^}w. (11.15) 

N N IN 1 1 


Now 


K -K 
N 


' (Ki„-Ki)-KK2^-K2) 


SO that the inner term in (11.15) becomes 

(K^ -K ) (H+K)“^+(K -K_) (H+K)~^-(K^ -K 1 (H+K 1 

IN 1 ^iN 2 IN 1 1 

= (K, -K. ){ (H+K)"^-(H+K^)”^}+(K- -K. ) (H+K. )“^ 
IN 1 1 2N 1 1 

= (K^jj-Kj^) (H+K^)"^(K^-K) (H+K)~^+(K2 j^-K2) (H+K) 


= { (K^^-K) (H+K^)'^(K^-K2)+(K2^-K2) } (H+K)~^. (11. 1G) 

Substituting (11.16) into (11.15), using (H+K) = ij) and taking norms gives (11.12)0 

Equation (11.12) shows explicitly how the convergence rate of to ij) depends 
on the splitting of K. If K 2 can be chosen so that it is very smooth, then on the 
basis of (11.12) we would expect rapid convergence of to ijj. (This probably justi- 
fies its success in three dimensional free air calculations.) At present the only fea- 
sible splitting when walls are present seems to be the trivial one, K^=0, K 2 =K. In 

this case ip is determined by 
o 

HiJ; = w (11.17) 

o 

which can generally be obtained analytically via the Sohngen inversion formula [30] . 
The residual pressure then satisfies 

(H+K)i[)^ = -Kh”^w, 

R 

and is obtained from 

(H+K )iii = w. 

N N 
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The error estimate 11.12 novj simplifies to 

il ^ 1 1! (H+y"^llK<i»-K^1'||^ . (11.18) 

Using (11.18) we can directly compare the rate of convergence of to that of 
For the rate is proportional to ||H)J>-L^(Hi|>) 1|^ while for it is propor- 
tional to [|Ki}<-L^(K!l;)j|^. Since Kt|< is continuous should converge faster than 

i|i , but for discontinuous downv:ashes of the type associated with flaps v;e anti- 
^ 2 

cipate a rate no better than 0(1/N ) since Kij; v?ill in general not be differenti- 
able due to the logarithmic singularity in K^. 

Thus Landahl's method can be expected to give some improvement, but at the 
expense of having to determine the functions and Generally KiJ;^ v/ill have 

to be calculated to high accuracy numerically. This requires the evaluation of the 
Jcernel K{x-^) at the points v;here {x.}^ are the collocation points and 

{gj}^ are quadratxire nodes. Using logarithmic Gaussian quadrature one should be 
able to evaluate these integrals fairly efficiently. Since for Bland's equation 
the bulk of the computing time goes into the calculation of the kernel, there should 
be an N v/here the amount of arithmetic needed to implement Landahl's method is less 
for a given accuracy than using collocation directly on (11.1). This tradeoff point 
will have to be determined experimentally and will be a topic of future investiga- 
tion. 

11.1.2 Landahl ll 

In order to achieve a possible analytic simplification in the evaluation of 

lb via (11.17) we consider the following modification of the previous procedure, 
o 

Assume that w can be decomposed as 

= w + w (11.19) 

s c 

v;here w^ is the singular part of w and w^ is the continuous part. In this case 

iJ; = U'\f (11.20) 

o s 

and the residual dovmv/ash is given by 

V7_ = w -KH (11.21) 

R C 

The residual pressure is obtained from 

(H+K)ij;„ = w -KiJ; . 

R c o 


( 11 . 22 ) 
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N 

Eq. (11.22) is solved by collocation to give an approximation to The ap- 

proximation to (j; is taken to be 

Using arguments similar to those above we arrive at the following theorem. 

Theorem 11.3. Let be given by (11.23). Then converges to jp and 

II 'J'-’I'nIIyIII ^11 - (11.24) 

Proof. See Ref. [58] for details. Note that from (11.24) the convergence rate 

of Landahl II should be the same as that of Landahl 1, but that ili should be easier 

o 

to evaluate 0 
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11.2 I teration 

A common procedure in niimerical analysis for solving an equation is to use 
iteration. For linear equations Picard iteration is a well-knovm method. Appli- 
cation of this to (11.1) determines a sequence of approximations to by the fol- 
lov;ing scheme. Let be an initial approximation to ij) and recursively define 

t!i , n>i by 
n — 


+Ki/) , = w. (11.25) 

n ii“JL 

since ||h 1( = 1 a sufficient condition for to converge to is that ||k|[ <1. 

The error is given by 


II 'f'-'l’nllY:! II ^""11 II 

Prom (11.26) one sees that the closer is to rp the more rapidly converges to 
i|». If is chosen as a collocation approximation to ip then on the basis of the 
above argument it seems reasonable that iterating on it vrould be a useful method for 
accelerating convergence. The basic difficulty that v;e encounter is the apparent 
problem of calculating more than one iterate, so that we restrict ourselves to con- 
sidering the effect of a single iteration. 

Before proceeding, v;e point out that our attention was dravm to this method 
by the v;ork of Sloan, Burn and Datyner [59] , and Sloan [60] , on the use of iteration 
for Galerkin's method for integral equations of the second kind. It was apparently 
first demonstrated in [59] that iteration on a Galerkin approximation gave super- 
linear convergence. In addition to Sloan et al the method has also been used by 
Phillips in conjunction v/ith collocation for equations of the second kind [61] . 

Since vje vjill discuss both collocation and Galerkin's method in this section a slight 
change in notation is introduced, vjill denote the collocation approximation to 

il) and tlij! V7ill denote the Galerkin approximation [5] . 

N 

— c 

Let tl) be the first Picard iterate of ib„ defined by 


Theorem 11.4. converges vmiformly to the solution i{) of 11.1. 
Proof. From (11.27) it is seen that 


(11.27) 


(11.28) 
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It now follows from the Sohngen inversion formula that H is an integral operator 
with Hilbert-Schmidt kernel [5] . Let K denote the kernel of H Then 

f/F 


Similarly 


= w{x)- 


-1 




This gives 


<i){x) = w(x)- 




-1 


K(x,o^(0<ac. 


llj(x)-li)^(x) = 


-I 


/ k(x, 5) [v(?)-’J.^(C)ld5- 


By the Cauchy-Schwarz inequality 

max (X) |<c|| 

-l£x;<l ^ ^ 

But II ’^"'I'^ll 0 so that Vj^(x) converges uniformly to i()(x) D 

From the uniform convergence of i()^ to tp we anticipate that the sequence 

should be better behaved than # and possibly free of the Gibbs phenomena present 

2 ” 

in L convergence. One also expects a more rapid rate of convergence than that for 
collocation alone. For equations of the second kind the results of Phillips appear 
to indicate this [61]. However we have not been able to prove this at present for 
Bland ' s equation • 

Using Theorem (4.5) of [5] as a guide we expect that for even moderate N that 
the collocation and Galerkin solutions should be close, and it is for Galerkin's 

method that accelerated convergence can be established. 

G G 

Let W , be the Galerkin approximation to t|) ( 5 ] , and observe that tp , satisfies 
N N 


G G 

HiJj +11 Kip = IT w 
N N N N 


(11.30) 


v;here is the operator of orthogonal projection onto Span 

first Picard iterate of Then <p^ solves 

N N 


H'J'^K.P^ = w. 

Applying to both sides of (11.31) gives 


(11.31) 


S ,G 

IT Hip +TI K Ip,. = IT W. 
N N N N N N 


(11.32) 
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Since ijj satisfies (11.30) we get 
N 

Using (11.33) in (11.31) shows that ij;® solves 


H^f+K(H V,H)iJ;® = w. (11.34) 

N N N 

Let Q„ = H ^TT, H and observe that by the unitarity of H that Q„ is the operator of 
N N N 

orthogonal projection onto Span This gives as the solution to 

3 

We now consider the convergence of ijjjj to \}i. First, note that by the same argu- 
ment that was used in Theorem 11.4 that iJj® converges uniformly to i|). However, here 
we can establish a superlinear rate of convergence. Our result is a generalization 
of Sloan's [59] , [60] for integral equations of the second kind, 
s 2 

Theorem 11.5 ih converges in L to tii and 
N a 


II II(H+KQjj)"^|I II K-KQ^II II if--Qj^i(-|lY. (11.35) 

Proof. First observe that ||k-KQ |1 0 [5]. From this it follows that (H+KQ„) ^ 

exists for N sufficiently large.. Thus 


This gives 


= (H+K) ^w- (H+Kp,j“^w 
N "N 


= (H+KQ^ )“^(KQ -K)t|;. 
N N 


(11.36) 


Now 

so that 

(KQ^-K) ^ = (KQj,-K) (KQj^-K) { I-Q^) 4> . 

Since Qy = 




(11.37) 
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using (11.37) in (11.36) gives 

= (H+KQ^)”^(KQ^-K) (l|j-Q^l()) (11.38) 

Taking norms on both sides of (11.38) gives the resultQ 

By a similar argument to that above it is easily shown that 

II II • (11.39) 

Comparing (11.35) and (11.39) shows that the rate of convergence of to ip 
is enhanced by the factor ||k-KQ^|| over that of i/)^. For steady problems the log- 
arithmic terms in K will be absent so that || K-KQ^|| should converge to 0 rapidly 

Q 

thus providing a more rapid rate of convergence than For unsteady problems 

it follows from Theorem (4.5) of [5] that 

[[ K-KQ^II = 0(1//n) , (11.40) 

so that for discontinuous downwashes it follows that 

II V'-^®||^ = 0(1/N), (11.41) 

a rate equal to that observed for collocation. As indicated in section 10, this 

3/2 

may very well be pessimistic, however a rate no better than 0(1/N ) is anticipated. 

This at present appears to be too slow a rate to achieve the desired accuracy of 
0(10 ^) error with N<24. 
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11.3 Reverse flow theorems 


A method of long standing in the aerodynamics literature for solving flow 
problems has been the use of reverse flow theorems [55] , [56] . Originally intro- 
duced as a technique for decreasing the computation time in calculating lift and 
pitching moment it also appears to have some promise for increasing the rate of 
convergence of integrated quantities for problems with flaps. 

It was shown in [5] that such quantities as lift, pitching moment and gener- 
alized aerodynamic forces could be calculated in terms of inner products <f, 
where f is the function representing the particular moment of interest. The basic 
idea of a reverse flow theorem consists of evaluating <f,ip>^ in terms of the solu- 
tion to an appropriate adjoint problem. 

To be more precise let H* and K* be the Hilbert space adjoints of H and K [5] 

2 

and consider calculating the inner product <f , where fEL^- Let ip* be the unique 
solution to 


(H*+K*)1(J* = f. (11.42) 

(That \p* exists and is unique follows from the fact that H+K is one to one and the 
Fredholm alternative. ) 

Theorem 11.6. <f,i|;> = <\p*,w> (11.43) 

Y o. 

Proof. By the definition of adjoint we find that 

<f,Tl)> = < (H*+K*) \!j> = <\l)* , (H+K) 4i> 

y Ct 

= <ih*,W> □ 

^ a 

From the theorem we observe that if f is "smoother" than w, then solving for 
ip* should present a better behaved problem than solving for ip directly if one only 
wants to calculate <f,rp>^ rather than i|; itself. 

Theorem (11.6) suggests the following scheme for approximating <f,ifi>_^. Let 
t);* be a collocation approximation to ij)* and take as an approximation 
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* ic 

to . It xs easily shown that converges to ijj so that converges 
to <w,if'*>^ = Thus the method is well defined. Using the Cauchy-Schwarz 
inequality it follows [48] that 


|<f,ijj>^-<w,iJ;*>^|<|| w|| ^{|| (K*-K^)i() 


f-f.. 


N" a 


}. 


(11.44) 


If f is a polynomial, then for N sufficiently large f = f^ and (11.44) becomes 




[ (K*- 


K*),|.^ 


] . 


(11.45) 


The convergence given by (11.45) is then analogous to that achieved using a 
"degenerate" kernel rather than a projection method [51] . 

Although (11.45) indicates good convergence for a "polynomial" moment the 
method discussed above will generally be ineffective if f represents a discontin- 
uous moment such as the hinge moment, because in this case the adjoint problem (11.42) 
is of the same type as (H+K)ijj=w. Since for some problems the reverse flow Theorem 
11.6 appears to have improved convergence properties we feel that its use should 
be further investigated. 


11.4 Other Methods 

There are other procedures that have recently been proposed for the solution 
of problems with flaps. Among these are the semi-analytic methods of Williams [62] 
and [63] , and several projection methods proposed by Milne [56] . These techniques 
have yet to be investigated and would require one to develop algorithms different 
from that employed in TWODI. Although we have ruled out none of the above as 
future candidates, our analysis indicates that none of these would provide suffi- 
ciently increased accuracy without making major alterations in our existing code. 
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11.5 Extrapolation 

As the discussion above pointed out it is presently felt that incorporating 
any of the above methods of convergence acceleration would not provide an acceptable 
degree of accuracy without making substantial changes in TWODI . For this reason 
we have begun investigating methods which can make the maximum use of the output 
of our program. As we presently have the capability of solving many problems at 
once, and with considerable efficiency, a reasonable approach is to try to combine 
these solutions effectively for enhanced accuracy. A standard procedure for doing 
this in numerical analysis is the method of extrapolation [50] . Here various solu- 
tions computed for different numbers of basis elements are combined either linearly 
or non-linearly to produce solutions which have higher order accuracy than any of 
the original ones. Such methods are commonplace in codes for solving differential 
equations [37] , [64] but seem not to have attracted much attention for the solution 
of integral equations [50] , [65] . As we shall see, their effective use requires 
one to have available known asymptotic expansions for the error in the numerical 
approximation as a function of the number of basis elements. Such expansions seem 
to be lacking in general for collocation methods and in particular for the solution 
of Bland's integral equation. However, using the error estimate (10.7) as a guide, 
along with the examination of results for the free air case we show how to obtain 
the proper form of the error and thus to extrapolate correctly. 

Although many forms of extrapolation are possible, we have chosen the analogue 
of Richardson extrapolation [50] because one can easily obtain the correct order 
of convergence. For other methods, such as Aitken extrapolation, this is not easily 
done [37] . V7e begin by presenting the basic theory of this method. 

Let a be a complex number and let be e sequence converging to a. Suppose 
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it is known that the error e = a-a has an expansion of the form 

n n 

e = a^ /n'^^+a /n'^2+. . . (11.46) 

n 1 2 

where ai< 02 < 03 * • • The sequence is rtot necessarily convergent but is assumed 

to have the property that 

m 

lim n^m(e - J a,/n°:) = 0. (11.47) 

II • -I 3 
n-»^ j=l 

Such an expansion is called asymptotic. 

For simplicity we restrict ourselves to the case where all the ' s are integers 

and in particular assume that a.=j. Thus e has the form 

i n 

e = a /n+a /n^+a /n^+ — (11.48) 

n 1 2 >3 

The basic idea of Richardson extrapolation is to compute subsequences 

and then to combine these in such a way to form new approximations which converge 

more rapidly to a than {a }. For example, consider -a . Then 

n n 2n n 


a-Z^ = [2{a-a^/ n-a /4n^ . . . ) - ( -a^ /n-a /n^ . . . ) ] = -a /2n^+. 
n 1 2 12 2 


(11.49) 


Thus {Z^} is accurate to order 1/n . {z^} is referred to as a sequence of first 

n n 

order extrapolations of 

For each j it is possible to calculate a sequence {z^} which has the property 


that 


„j 1 , j+1 2 ,3+2 

a-Z = a , /n +a . /n 

n J 3 


(11. 50) 


Such a sequence will be called a sequence of order extrapolations of {a } . 

... 2 . ^ 

For example, one possibility for Z^ is 


Z = n /3-2a_ +8a. /3, 
n n 2n 4n 


(11.51) 


and for Z' 


n 


Z = - a /21+2a^ -8a, /3+64a„ /21. 
n n 2n 4n 8n 


(11.52) 


I I III 
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Higher order extrapolation formulas can also be obtained, 
a derivation of (11.52). 

From (11.48) it follov/s that 

2 3 

e = a^/n+a_/n +a /n + — , 
n X ^ 3 

2 3 

= a^/2n+a2/4n +a^/8n + 

2 3 

e, = a, /4n+a_/16n +a /8n +. . . , 
4n 1 2' 2 


2 3 

E = a /8n+a /64n +a /512n +. , 
8n X 2 3 


(ct,g,Y) are now determined so that 


1 4 

E + ae_ +Y£n = 

n 2n 4n 8n 4 


From (11.53) this requires that 

a/2+8/4+Y/8+l = 0, 


a/4+S/16+Y/64+l = 0, 


a/8+6/64+Y/512+l = 0, 


We illustrate these with 


(11.53^) 

(11.532) 

(11. 533) 
( 11 - 53 ^) 

(11.54) 

(11.55^) 

( 11 . 552 ) 

(11.553) 


Solving (11.55) gives 


a = -14,6 = 56, Y = -64. (11.56) 

Substituting into (11.54) gives 

(a-a )-14(a-a„ )+56(a-a. )-64(a-an ) 
n 2n 4n on 

= a- (-a /21+2a„ /3-8a . /3+64a„ /21) = a^/21n^+. . . 
n 2n 4n 8n 4 


Letting 

= -a /21+2ct^ /3-8a^ /3+64n /21 

n n 2n 4n 8n 

shov;s that (11.50) is satisfied for j=3. 

In practice it is most common to use the subsequences (a 2 k obtain 

n 

the sequence {Z"^}. Solving the following set of linear equations allovis one to 
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obtain the appropriate coefficients for jth order extrapolation. 


1/2 1/4 ... 

1/4 1/16 

f 

CN 

H H 


'^l' 


r ^ 
1 

1 

1/2^ 1/2^^ . . . 

n2 

1/2^ 

a 




1 


(11.57) 


In this case we have 


E +B,e +BtE. + — +6.E_j = a.,T, 

n 1 2n 2 4n 3 2 n 3 +I 


/n +. . . 


(11.58) 


Since the coefficient matrix in (11.57) is a Vandermonde matrix [5] , equation 

(11.57) is uniquely solvable and the sequence {2^} satisfying (11.58) is well defined 

To determine the applicability of these formulas to the solution of problems 

with flaps the results of Tables 3-6 were used. Here we observe the curious fact 

that for a midchord flap the subsequence {a_ (a = lift or norm) appears to sat- 

ri“*i n 0 

isfy (11.48) and for a three-quarter chord flap the subsequence appears to 

have the same property. Thus for a mid-chord flap it seems to be appropriate to 
extrapolate on {a-,a. ,ao/a. and for the three-quarter chord flap we use 
This oscillatory behavior, with the period of oscillation apparently depending on 
the location of the hinge point was unanticipated based on our experience with the 
rapid and monotone convergence of TWODI-I for smooth downwashes. At present it 
appears that the period of oscillation is independent of the kernel K so that a 
suitable subsequence for extrapolation can be obtained by solving a series of steady 
free air problems . This can be done rapidly with our existing program. 

1 2 

To see the effect of using successively higher order extrapolations Z , Z , 

3 1 ° 2^ 

and Z were calculated for the lift and norm for the mid-chord flap and Z , z 
2 6 3 

were obtained for the three quarter chord flap. The results are listed in Tables 

7-9 below. 


Table 7. Extrapolation of and noimi for a midchord flap (M=0,k=0, 



Norm 

1.60337 

^^6 

5.14159 

% error 

% error 

Exact 

(norm) 


^8 

1.60241 

5.13466 

-.0599 

1348 

2 ^ 

1.60327 

5.14069 

-.0062 

-.0175 

Z^ 

2 

1.60335 

5.14139 

-.0012 

-.0039 
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Table 8. Extrapolation of for an oscillating flap hinged at the 50% chord 


(M=0, k=.l. n=o») 


real 

imag 

magn 

phase 

% error (magn) 

Exact 

4.33227 

-0.43177 

4.35373 

-5.69 



4.32647 

-0.43122 

4.34791 

-5.69 

-0.1337 


4.33153 

-0.43171 

4.35299 

-5.69 

-0.0171 


4,33210 

-0,43176 

4.35359 

-5.69 

-0.0033 


Table 9. Extrapolation of 

for an oscillating flap hinged 
(M=0, k=.l, n=“) 

at the 75% chord 


real 

imag 

magn 

phase 

% error (magn) 

Exact 

3.20444 

-0.48982 

3.24166 

-8.69 


6 

3.19235 

-0.48792 

3.22942 

-8.69 

-0.3776 

2 

Z 

3 

3.20191 

-0.48944 

3.23910 

CD 

VD 

-0.0790 


Table 

10. Effect of 

usina incorrect 

subseouences for extrabolation 

real 

midchord flap 

imag 

magn 

phase 

% error (magn) 

Z^ 

6 

4.32230 

-0.43083 

4,34372 

-5.69 

-0.2300 

2^ 

4.57927 

-0.45301 

4.60163 

-5.65 

+5.6939 

three-quarter chord 

flap 




Z^ 

8 

3.51114 

-0.53452 

3.55159 

-8.66 

9.5608 

Z^ 

4 

3,77727 

-0.57288 

3.82352 

-8.92 

17.9494 

Z^ 

2 

4.00144 

-0.60507 

4.04693 

-8.60 

24.8413 
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The results are quite interesting showing that 4-5 decimal accuracy is obtain- 
ed using at most 16 basis elements. Thus it appears that an error of 0(10 can be 
obtained using the present code, followed by the minimal cunount of manual arithmetic 
required to do the necessary extrapolations. 

As an independent check on the validity of the above procedures we have extra- 
polated some results presented by Milne [56] in his development of finite element 
methods for solving (3.8). Table 11 taken from [56] gives values of lift, pitch- 
ing moment and hinge moment for approximate solutions to (3.8) using basis elements 
which were piecewise linear. 

2 

Again it is easy to show that the errors £ satisfy e = a,/n+0(l/n ) so that 

n n 1 

extrapolation is feasible- Second order extrapolations are given in Table 11. It 
is interesting to note that the "exact" values are obtained using 48 basis elements 
with extrapolation, whereas the raw value for lift and hinge moments are correct 
to only 2 significant figures using 60 basis elements. 


Table 11. Convergence of section coefficients for 50% flap using finite elements 

(M=0 . k=0 . n=”) 


No. of basis 



Pitching 

Hinge 

elements 

Lift 


Moment 

Moment 

Exact 

-3460 


.1365 

.0780 

12 

.3410 


.1353 

.0756 

24 

.3435 


.1359 

.0767 

36 

.3443 


.1361 

.0771 

48 

.3447 


.1362 

.0775 

60 

.3450 


.1362 

.0780 


.3460 


. 1365 

.0780 

12 





On the basis of these observations 

we 

feel confident that the use of extrapola 

tion will result in accurate 

solutions 

for 

flow problems 

with flaps. The main ob- 

Stacie remaining to complete 

automation 

seems to be the 

necessity of performing the 


preliminary analysis of a sequence of free air problems in order to determine the 
appropriate subsequence for extrapolation. In this regard it is interesting to note 
that while Bland's collocation method produces oscillatory convergence, Milne's 
finite element procedure does not. Monotone convergence is, as we have seen more 
desirable and further investigation of the methods studied in this section for this 
behavior should be done. 
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§12., Efficie nj^ integration of singularities in the kernel 

The basic accuracy of the TWODI program for smooth downwashes is about 6 decimals 
in steady flow and 3 decimals in unsteady flow using 5-10 basis elements. We will 
now show how this can be improved to 6 decimals in both steady and unsteady flow with 
little increase in the amount of arithmetic. 


12.1 General procedure for logarithmic singularities 

An efficient solution of the airfoil equation (3.16) requires that the integral 
for the collocation matrix 


C 

mn 


il 


K(x - 
1+^ m 




.)Tn(C) 


d? 


( 12 . 1 ) 


be computed efficiently for all values of the parameters etc. In general, the 

kernel contains a dominant Cauchy singularity, a weaker logarithmic singularity, and 
is otherwise analytic. It can be shown that 


K(x,k,M,n, 


.) = 


ik -ikx„ II -ikx „ , , „ 

e Fj (x,k,M) log|x|-e F 2 (x,k,M,n». 


82 


.) , 


( 12 . 2 ) 


where F-]^ and F 2 are analytic. We will identify these functions below. 

The procedure originally used by Bland [ 8 ] , [ 9 ] to compute is the one 

presently used in TWODI and is described in section 6. The Cauchy singularity is 
integrated in closed form using the Akheizer-Bland transform (5.3) and an additional 
closed form integration is obtained by using the logarithmic transform (5.24). The 
remaining continuous part of the kernel multiplied by the appropriate basis function 
is then integrated approximately by Jacobi-Gaussian quadrature (5.30). Thus the 
only source of error (other than roundoff) in computing Cmn presently is quadrature 
error from integrating the continuous part of the kernel. 

It is well known [37] , [38] that high precision Gaussian quadrature formulas 

work best with analytic integrands. Splitting the kernel into the Cauchy, logarith- 
mic and continuous parts as described in section 6, 


1 


^ log|xi + (x,k,M,g, . . . ) , 


K(x,k,M,T|, . . . ) 


X 


(12.3) 
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v/e see by (12.2) that the continuous part of the kernel is given by 

K|^(x,k,M,n, . • . ) = (x,k,M)) log|x] - e (x,k,M,ri, . . . ) (12.4) 

aind is clearly not analytic because it contains nondif ferentiable terms of the 
form 

X log I X I , log | x | , . . . 

To assess the errors resulting from integrating without properly accounting 
for these logarithmic terms, we computed 

1 

x^^log X dx 

^ 0 

using Legendre-Gaussian quadrature for various n and for various numbers NQ of 
quadrature points. The results, shown in Table 12, indicate that for n=l, 16 point 
quadrature is at best 5 decimal accurate. 


Table 12. Error in j 

^ n 

X logxdx 

1 0 

using NQ point Legendre-Gaussian quadrature 

% error 

n = 0 

n = 1 

n = 2 

NQ = 1 

30.685 

-38.629 

-55.958 

NQ = 2 

10.412 

-3.1413 

2.3139 

NQ = 4 

3.1464 

-.25967 

.04269 

NQ = 8 

.87610 

-.01956 

.00083 

•zi 

JO 

II 

H 

O 

.23207 

-.00136 

.00001 


Since the number of quadrature points used in TWODI equals the number of pres- 
sure basis functions, it is not reasonable to expect to obtain accurate calculations 
for unsteady flow using a small number of basis functions. For example, using 8 ba- 
sis functions we might expect at most 3 or 4 decimal accuracy. This has been borne 
out in practice; using NP = 8, the TWODI-I program was 3 decimal accurate for un- 
steady flow [5, p.81]. 

In order to correctly integrate the logarithmic kernel singularity in the pres- 
ence of leading and trailing edge pressure singularities, it is necessary to separate 
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the interval [-1,1] into four stibintervals [-l,a] , [a,x] , [x,b] and [b,l] where 

-l<a<x<b<l. (12.5) 

Making the necessary transformations gives six different kinds of integrals, all 
of which can be evaluated quite accurately (12 decimals) using tables stored in 
TWODlJ 


n 


-1 


l+c 


log |x-?| f(x,E)d5 


= >^l+a 1 

TT J 0 /S' 


[/l-? log (x-C) f (x, C) ] du 


x~a 

TT 


0 u f(x,?)ldu 


(x-a) log (x-a) 


n 


/i-c 


[ / TTF f (x,5) ]du 

0 1+5 


(b-x) log (b-x) 

IT 


+ ^ ^ I [ /— p f (x, 5) ]du 

1 0 '' 1+5 


n f/i 


-5 


b-x 

TT 


0 [/f^f(x,?)]du 


'l-h 


n 


1 r 1-5 


;= [ 1 =^ log (^-x) f (x, O ]du 
0 /u /1+5 


5 = -l+(l+a)u 


5 = x-(x-a)u 


? = x-(x-a)u 


C = x+(b-x)u 


5 = x+(b-x)u 


? = l-(l-b)u. 


(12.6) 


Each of the six integrals appearing in the right hand side of (12.6) possesses an 
integrand which is the product of an analytic function shown in square brackets, 
multiplied possibly by a function representing an inverse square root or logarith- 
mic singularity. For obvious reasons, we call these six integrals the leading edge 
inverse square root, upstream logarithm, upstream continuous, downstream continuous, 
downstream logarithm, and trailing edge square root parts, respectively. 
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12.2 The free air incompressible case 

This is the simplest case with which we can test the merits of integrating 
the logarithmic singularity in the kernel using high precision Gaussian quadrature. 
We will show that a thousand fold increase of accuracy results, with only five ba- 
sis functions being required to produce six decimal accuracy. Previously, eight 
basis functions produced three place accuracy. Since accuracy can be traded for 
reduced computer time, the results of this subsection constitute an important step 
toward TWODI handling its own accuracy control in the future. 

In this case the kernel may be written as [5,p.7] 

K(x,k,0,co,...) = i - ike (Ci(k X ) + i Si(kx) + ^) (12.7) 

X 2 


where Ci and Si are the cosine and sine integrals. They are given [66, p. 232] by 


Ci(z) = Y + log z + ^ 


(- 1 ) 


n 2n 


n=l 


(2n) (2n) I 


( 12 . 8 ) 


Si (z) 


1 


n=0 


, . n 2n+l 

(-1) z , 

(2n+l) (2n+l ) I 


(12.9) 


where y =-57721566490153... is Euler's constant. Combining the aibove, we obtain 


K(x,k,0, »,...) = ^ -ike log[x| - ike“^^^(log k + y + ^ + I (12.10) 

n=l ' 


Thus we can identify the analytic functions in this case as 

F^(x,k,0) = 1, 


( 12 . 11 ) 


(x,k,0,«. 


. . ) = ik (log k + y + ^ 

m=l 


(ikx)^ 
(n) (nl)^ 


( 12 . 12 ) 


Tlie effect of using equations (12.6) and (12.10) to integrate the complete log- 
arithmic singularity, leaving an analytic integrand to be done by Jacobi-Gaussian 
quadrature is shown in Table 13 for the case of a plunging airfoil. The exact solu- 
tion was obtained from the Kussner-Schwarz comparison [5,p.83]. Using the earlier 
method based on the logarithmic transform with anxlog|x| singularity in the contin- 
uous part of the kernel, 8 basis functions produced only 3 place accuracy. Table 13 
shows that this is true both of TWODI-I using Bland's kernel with very large height 
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to chord ratios (and large computer times) , and of TWODI-III which calculates the 
special kernel (12.10) which is more efficient; this verifies that using Bland's 
kernel with n = 300 is an accurate representation of (12.10). Thus the rather low 
3 decimal accuracy obtained by the original method is due to the manner of integrat- 
ing the remaining part, K , of the kernel. Finally, Table 13 shows that using the 

c 

analytic form of the remaining part of the kernel produces six decimal accuracy with 
only 4 or 5 basis functions. 

The evidence presented by these data strongly indicate that the mos-t efficient 
computing strategy dictates a switchover to integrating the singularities in the 
kernel separately. 


Table 13 . ^ Method of integrating singu!^arit ies in^ the kernel 


Method 

type 

NP 

Cli “ plunge 

mode (k = 1) 

Exact 

- 

- 

-2.51156+3.389371 

= 

(4.21850,126.54°) 

TWODi-I 

X log 1 X 1 ^ 

8 

-2.50990+3.387581 

= 

(4.21608,126.54°) 

TWODI-III 

X log|xp 

8 

-2.50990+3.387571 

= 

(4.21607,126.54°) 

TWODI-III 

analytic^ 

1 

-1.59997+1.584141 

= 

(2.25153,135.28°) 

TWODI-III 

analytic^ 

2 

-2.55260+3.426141 

= 

(4.27249,126.69°) 

TWODI-III 

analytic 2 

3 

-2.51089+3.389001 

= 

(4.21781,126.53°) 

TWODI-III 

analytic 2 

4 

-2.51156+3.389371 

= 

(4.21851,126.54°) 

TWODI-III 

analytic 2 

5 

-2.51156+3.389371 

= 

(4.21850,126.54°) 


^Using Bland's kernel with n = 300. Approx. 20 min CPU. 


^Using eq. (12.10) for kernel. A few sec. CPU. 
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12.3 The free air compressible case 

Identification of the singularities in the Possio kernel is certainly impor- 
tant to an efficient solution of the free air compressible case, and also may prove 
useful in efficiently computing more general kernels for unsteady flow in wind tun- 
nels with slotted and porous walls. 

The Possio kernel may be expressed [5,p.9] as 


K(.,k,M, ^ ,(a^ (iMsgnW 


2i3 , 1+3 ^ ^ ikX „ (2) .kMlXl ^ 

+ — log — + rk j ^ e ^ Ho (~^) ) dA , 


(12.13) 


where the Hankel functions of the second kind [66, p. 360] are given by 


Clearly (12.13) is not in a form suitable for numerical computation. Substituting 
(12.14) gives 


K(x,k,M,“, . . . ) = e 


-ikx 


ikx 

'23^ 


TTkMi 

23^ 


sgn ( X ) 


(Jl( 


k M I X I 
32 


-) - i 




k M|xl 

32 


■)) 


+ 


TTk 

232 


(J ( 

o 


kwIxL 

32 ’ 


iY (■ 
o 


kM 


32 


-))] - 


3 


log 


1+3 

M 


ikxu 

TTk2xi 32 

232 Jo ® 


J (Mi^) 

O 


iY (Ml|lil))du}. 
o 32 


(12.15) 


In order to identify the singularities in (12.15), we first identify the speci- 
fic singularities, as well as analytic functions, which appear in the Bessel func- 
tions. Clearly 


J,(z) 




m! (m+5,) I 


(12.16) 


is analytic. Let 

1 1 

ii)(D = - Y, ijj(n+l) = rlj(n) + - = - y+I - if n>^2 (12.17) 

m=l 

denote values of the psi or digamma function [ 66 ] for integer arguments. Then 
it follows that 
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Y (z) = — (Y+log^) J (z) + — G (z) , 
O TT 2 O IT O 


Y, (z) = - — + — log -f- J (z)- — G (z) , 
1 itztt 21 irl 


where 


1 f)" 1 1 

^^-^2 3^— 


(in 


( 2 !) 


(3!) 


G (z) = J I (iMn+l)+ijj (n+2) ) 

^ ^ n=0 n! 


z 2 n 

<-| ^ 


(n+1) ! 


Obviously G^ and G^ are analytic. We note that they may be computed using 
or, using standard Bessel function subroutines according to 

G (z) (z)-(Y+log^) J (z) , 

O 2 o 2 o 

G-, (z) = ^ Y (z)-— + log^ J (z) . 

2 ^ z 2 1 

Combining the above and rearranging it can be shown that 
, 1 ik -ikx , I , kMx > , kMX ' , ikx 

K(x, k,M, ~ e log[x| [ (J^(-j2")- ^ ®"b2" 


- ikx 


i: 


ikxu , ,kMxu^ ^ , 
® '12-^ 


ikM^x 


-ikx, ,e b 2 
e t [- 


X 

kM 


•nk kjytx . ik , kMx 


ik , . kM , _ ,kMx, . irkMi ^ ,kMx, 

+ (Y+log —;pr^ J- (^7^) + J-, 


p2 


262 1 6" 
ikx 


kM , kM ^ , kMX ' kM , kMx , .. 

* 6? 5bT 62 h'-p-H'-e" 


ikxu 


ik , 1+6 . TTk xi _ 62' 

Jo ^ 


,kMxu. 2i „ ,kMxu. . , 


k2x , , . kM ^ 

-p- (Y+log 


1 


ikxu 

. 62 


0 

ikxu 


^ ,kMxu^ , 

Jq( g2 > 


k X 
62 


1 ft2 , ,kMxu, , 1 
log — e P J ( ) dui , 


o' 62 


(12.18^) 

(12. 18^) 

(12.19^) 

(12.192) 

(12.19) 

( 12 . 20 ^) 

(I 2 . 2 O 2 ) 


(12.21) 


0 
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which in view of ( 12 . 2 ) gives 


ikx 

F^(x,k,M) = - iM J^(^)]e - ikx 


i: 


ikxu 

32 kMxu, 

e J (■ o ' o ) du 
o 


( 12 . 22 ) 


and 


(x,k,M,«>, . . . ) = 


- i3cM^x 
e P 2 -1 


TTk , kMx , , ik , kMx , 

2^ 


. ik , kM , , ,kMx, , TTkMi ^ ,kMx, 

+ -p- (Y+log ^) J^(^) + 


262 o 6" 


232 


kM , kM ,kMx kM ,kMx, ^ o 


ikx 


ikxu 


^ T IT 


Trk^xi 

fl - 

62 

(J ( 
0 

kMxu. 
B 2 > 

2 i 

262 

e 

0 

TT 

ikxu 





kM . ^ 

^ Jo ^ 

62 

J 

0 

,kMxu.^ 
'■ 62 ‘ 

du 



. ,kMxu 
62 


■) ) du 


, f 1 ikxu 

k X , 1 ~~St~ t /kMxu. , 

log - e J (-^^) du. 

62 J 0 u o 62 


(12.23) 


Equations (12.22) and (12.23) display the desired functions F^ and F^/ but 

they are not yet in final form and (12.23) is still indeterminate for M=0. F^ 

can be put into final form by first noting that 

J- = - (12.24) 

o 1 


and integrating by parts to yield 
ikxu 

ri 

ikx 


62 


.kMxu. , 


„o , .kMx, 

= B^(J (-^)e 


ikx 

62 


-1) +kMx 


ikxu 
. B 2 


, ,kMxu, 


(12.25) 


Substituting (12.25) into (12.22) and simplifying, we finally obtain 


2 

F., (x,k,M) = 1+M (J (-^)e 
1 O p2i- 


ikx 

62 


/kMx. 

1) - iMJ (-^)e 
1 p 


ikx 

62 


kMx 


ikxu 

^ e 
0 


.kMxu. J 


( 12 . 26 ) 


In this form (12.26) reduces by inspection to (12.11), is easy to compute, and we 
see that 

F^(x,k,M) = 1+M^0(kx) (12.27) 

is analytic as asserted earlier. 
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We now obtain a similar reduction for F^. Integrating the last term in (12.21) 
by parts gives 


k2x 1 


ikxu 

62 kMxu 

e du 


ikxu 


" Jo = ^o<— > + -p- J 

(We note that when M=0, (12.28) reduces to 

■ 1 ikxu 


n 


. 1 , ikxu , ^ ,kMxu, 

0 u ® 62 “ ■’^l ~62~^ (12.28) 


ik 


-1 


0 u 


du = ik 


n=l 


(ikx) 

(n) (n! ) ' 


(12.29) 


which is the last term in (12.12).) Substituting (12.28) into (12.23) gives 

ikxu 

•1 


(x,k,M,“>, . . . ) = ik 


-1 ,3cMxu, , 

J du 


ik^Mx 

B2 


1 
0 

ikM2x 


u o 

ikxu 


log i (e - 1 ) du 


62 


+ ik[- 


-1 . Itt .. ,kMx 


kMx. 


ikx 


2 62 +i2 


. 1 . 1 /. .g kM . ,kMx, ttM ,kMx. 

62 262^ o‘“6^^ 262 '^1^'^^ 


ikx 

kM , kMx , . iM , kMX' , 62 


.Ik 1+6 Trk^xi 

T “5“ * 


ikxu 

” e (cr (igH). lie te, d„ 

O 77 O $2 


+ ^ (y+log ^) 


n 


ikxu 


„ 62 ^ ,kMxu, , 

0 ^ 


(12.30) 
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Integrating by parts again, it follows that 

ikxu 


^ ,kMx, xkx 
o 6 ^ B 


ikxu ikx ikxu 

B2 ,kMxu, , -r /kMx, 32 , ,, 32 , ,kMxu, ^ , 

e J (— )du=3 J (-r^)e -kMx e J., ( -■ , )du. (12.31) 
Jo o 32 o 32 Jo 1 32 


Substituting (12.31) and the identity 


, kM . irk , ik 1+3 ■! /i , , ik , -1 „ 23^ , , 1+3 ' 

ik(Y+log — + — log — = xk(logk+Y+ ~ ^1+3^°^ ~ 


232 2 ' 3 

into (12.30) gives 


ikxu 


F 2 (x,k,M,“>, — ) = ik[logk+Y+ ~ + 


e -1 


M 


^ ,kMxu, , , 
J ( ) du] 

o 3^ 


232 ' 


(12.32) 


-ikM2x 

32 9 

,e -1 ,iir kM ' ,kMx. 

^ ^ ^ 


iM.iiT kM ' , ,kMx 

^(-+iog^) J^(-^) 


ikx 


1 _ , kMk ' , iM , kMx , 3^ 
g 2 32 


ik 2 > 

32 


ikxu 


■‘r.., kM . iTT' ^ ,kMxU',.„ ,kMxu, 32 

0 [M (Y+log^ + — ) (-^) +iG^ (-^) e du 


+ 


ik^Mx 

" 3 ^ 


1 

log — 
0 u 


, ,kMxu, , 


ikxu 

32 


- 1 ) 


du 


^ ik,M2 232 

* T'Tfe 




(12.33) 


which, in view of (12.29), reduces by inspection to (12.12) when M=0, and is clearly 
analytic with the property that 


F 2 (x,k,M,“, . . . ) 


F 2 (x,k, 0 , . . . ) 


+ 0 (kM log M) . 


(12.34) 
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One reason for the importance of (12.27) and (12.33) is the ease with which 
they can be computed niamerically compared with the original expression (12.13). 

In TWODI, Bland's kernel (3. 10) - (3. 11) is presently used for all cases except the 
incompressible free air case given by (12.10) and the steady porous wall case 
given by (4.13). However, using Bland's kernel for free air conditions is ineffi- 
cient because the series expressions in (3.11) are quite slowly convergent for 
large values of g. We have already seen in Table 13 the dramatic benefit of using 
the efficient expression (12.10) to compute the kernel for the special case of in- 
compressible flow in free air. We would expect a comparable reduction in computer 
time if (12.2), (12.27) and (12.33) were used to compute the kernel for compressible 

flow in free air. Furthermore, (12.2) represents a computationally efficient split- 
ting of the kernel into three parts 


K(x,k,M, rir • • ■ ) = — + K (x,k,M) + K- (x ,k,M, n , . • • ) , 

X 1 2 


(12. 35) 


where 

K^(x,k,M)= - -p- log|x| F^(x,k,M), (12.36) 

and where 


K2(x,k,M,n,...) = K2(x,k,M,», 


. . ) + AK^ (x,k,M,r|, 


) 


(12. 37) 


with 


(x,k,M,“, . . . ) 


-ikx „ , , V 

-e (x,k,M,“, . . . ) 


(12.38) 


and with AK 2 representing the interference kernel due to the wind tunnel walls. Since 
we now have computationally efficient expressions for all quantities except AK 2 and 
since the Fourier transforms of all kernels are easily obtained, it may be that the 
application of Laguerre-Gaussian quadrature to the expression 


AK 2 (x,k,M, n. 


) = 


/2tt 


^ g V A A 

5 [K(s,k,M,Ti, ... )-K(s,k, M, «>,...) Ids 


(12.39) 


is an efficient way to compute numerically the interference kernels for porous as 
well as slotted wall tunnels. Since the determination of the unsteady porous wall 
kernel remains an open problem, the possibilities of this approach will have to a- 
wait future research. 
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§13. Steady airloads for porous wall tunnels 

In this section we present numerical computations for steady airloads in ideal 
porous wall wind tunnels; i.e. , those governed by the viscous effect boundary con- 
dition (2.10). Agreement with previous results is demonstrated and new results 
are given in the form of compact three dimensional graphs showing airloads for all 
possible combinations of Mach number, height to chord ratio and viscous effect ven- 
tilation coefficient. These graphs are for a fixed downwash but otherwise include 
all possible flow conditions. They can be generalized to unsteady flow but we have 
not done so. This section is an extension of similar results by us [5, §§11-12] for 
fixed Mach numbers in ideal slotted wall tunnels; i.e., those governed by the mass 
effect boundary condition (2.9), and uses the parametric analysis presented above 
in section 4 to combine three independent variables (M,ri,v) into two (Bn,Bv). 

Table 14 and Figure 4 present values of lift coefficient and center of pressure 
vs. leakage angle The center of pressure is relatively insensitive to small 

amounts of ventilation and shifts more strongly forward as the open jet condition 
is approached. The trends in both lift and center of pressure using the ideal 
porous wall boundary condition differ appreciably from the trends obtained using 
the ideal slotted wall condition. This is evident upon comparing Figures 4 below 
with Figures 14 and 16 in [5] . Since the limiting case of open or closed walls are 
independent of the boundary condition selected, verification with previous results 
is obtained automatically. Only 5 pressure basis functions were required to obtain 
6 decimal accuracy. Also, the computations proceeded quickly since the kernel 
(4.13) is in closed form. 

Table 15 and Figures 5, 6 and 7 present lift coefficient, pitching moment co- 
efficient and center of pressure for the full range of Mach number, height to chord 
ratio^ and ventilation. The capability to do this so simply is a direct result of 
the parametric reductions 

(x,M,n,v) ^ (|^, Bn,6v) 

presented in section 4. The doubly infinite domain of the airload surfaces is made 
compact by the transformations 

^ = tan-l(^). 


^Tn the present context, we restrict our consideration to tunnels with acoustic 
height to chord ratios not less than 1; thus l£Bn£“>. This restriction is entirely 
one of context since the numerics are well behaved for narrower tunnels. 


-109- 


The lines 0-^=0 correspond to the free air condition and show constant airloads for 
all ventilation conditions; this is in keeping with the physical condition that air- 
loads not be affected by the tunnel walls when they are infinitely far away. The 
lines ?=0° correspond to a completely closed wall and indicate that both and 
increase as the walls become closer together; at the same time the center of pressure 


moves aft. The lines C=90° correspond to an open jet tunnel; Cx,^ and decrease 


as the walls become 

closer while 

the center 

of pressure 

shifts forward. Between these 

bounding 

lines , the 

airloads are 

shown to be 

continuous. 

Finally, by comparing the 

lift and 

moment surfaces, it is 

clear that the ventilation coefficients corresponding 

to zero 

lift interference differ 

from those 

corresponding to zero moment interference. 

Table 

14. Lift and center of 

pressure vs. 

ventilation 

for a Dorous wall tunnel 


M=.85 

n=7.5 

w (x) =1 

NP<5 (6 decimal accuracy) 



C 

Ct 

ija 

^CP 




0° 

12.235 

.25624 




15° 

11.460 

.25633 




30° 

10.739 

.25529 




45° 

10.051 

.25529 




60° 

9.4086 

. 24936 




75° 

8.8017 

.24416 




90° 

8.2270 

. 23723 



Table 15. 

Section coefficients vs 

. M, n and v 

for unit downwash 


3 Ct 


6 = tan~^ (-^) 

n 6n 



r. 

0° (free) 

15° 

30° 

45° 


0° 

6.28319 

6.46436 

7.06254 

8.29957 


30° 

6.28319 

5.62903 

5.31359 

5.28501 


60° 

6.28319 

4.89657 

3.97490 

3.34116 


90° 

6.28319 

4.24753 

2. 89934 

1.91357 




®n = 



C 

0° (free) 

15° 

30° 

45° 


0° 

. 000000 

-.022470 

-.094292 

-.233113 


30° 

. 000000 

-.016798 

-.076374 

-.197892 


60° 

.000000 

.001509 

-.011706 

-.060754 


90° 

.000000 

.030282 

.085436 

.132578 


^CP 



(-^) 

Bn 



X. 

0° (free) 

15° 

30° 

45° 


0 ° 

. 250000 

. 256952 

.276702 

.306175 


30° 

. 250000 

. 255968 

. 278747 

.324888 


60° 

.250000 

.249384 

.255890 

. 286367 


90° 

. 250000 

.235741 

.191065 

.111435 
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§14. Steady and unsteady interference calculation^ _for airfoils with flaps 

The nature of unsteady aerodynamic interference in wind tunnels is not well 
understood at the present time, especially when the walls are ventilated. This 
problem is made mathematically more difficult when flaps are present because of 
the attendant downwash discontinuity at the hinge and the resulting slow conver- 
gence of the numerical calculations. In this section, we present new results for 
steady and unsteady lift interference on airfoil with flaps at various combinations 
of Mach number, reduced frequency and wall ventilation condition. These results 
are based on the method of convergence acceleration by extrapolation to the limit 
as described in section 11 above. 

We take as our configuration of interest an airfoil of chord 180mm with a 45mm 
trailing edge flap mounted in a wind tunnel of height 550 mm. All calculations be- 
low were perfoiatied using second order extrapolation with NP=3,6 and 12. Based on 
our analysis in section 11, we expect these results to have mathematical errors of 
less than .1%. 

Table 16 and Figure 8(a) show the envelope vs. Mach number, along with the 

exact free air solution. The Mach niambers listed correspond to the Multhopp angles, 
M=sin 6, 0=0°, 18°,..., 72°. Figure 8(b) shows the interference ratio for Cl^ vs. M 
using the same data. In all cases the effect of wall interference is to increase 
the lift for a closed wall and to decrease it for an open jet. At very low speeds, 
these effects are roughly +3% and -36% respectively, increasing in magnitude as the 
Mach number increases. The increase is less and is more delayed with the closed wall 
and it is greater and more gradual with the open jet. At high subsonic speeds, es- 
pecially those above the critical Mach number, transonic nonlinearities will increas- 
ingly alter the interference effects from those predicted by the present linear theory. 
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Table 16. Steady lift interference vs. M for a flap at the three quarter chord 



n=3. 05556 


+ 

^ (exact) 


closed wall 

Open 

jet 

M 

Ct 

C 


p pfree 

free 


free 

.000000 

3.93081 

+.027274 

2.45693 

-.357908 

.309016 

4.14457 

+.030125 

2.53450 

-.370054 

.587785 

4.92542 

+.014371 

2.78953 

-.410216 

.809017 

7.00387 

+.075873 

3.30696 

-.492014 

.951057 

15.26128 

+.232474 

4.35129 

-.648597 


The phenomenon of acoustic resonance was orgininally discovered theoretical- 
ly by Runyan and Watkins [67] in 1953 and subsequently verfied experimentally by 
Runyan, Woolston and Rainey [36]. Later, Bland [8], [9] showed that for ideal 
slotted wall tunnels described by the boundary condition (2.9), acoustic reso- 
nance will occur at reduced frequencies given by 

MX 

’"n = n=l,2,... (14.1) 

where X^ is the nth positive eigenvalue of ( 3 .II 4 ). Theoretical calculations 
showing the effect of acoustic resonance on airloads over the full range of 
ventilation coefficients were recently given by Froinme and Golberg [5] . 

Table 17 and Figure 9 present values of vs. Mach number for open and 
closed tunnel walls at four values of reduced frequency, k=0,.l,.2 and .3. The 
calculations utilize Bland's kernel (3.10) for an ideal slotted wall based on the 
mass effect boundary condition (2.9). A condition of acoustic resonance between 
the airfoil and the wind tunnel walls is displayed at k=.249 for M=.9. A pre- 
cursor of acoustic resonance may be detected at M=.8 in that the magnitude of the 
lift coefficient for the closed wall condition has dropped below the value for 
the open jet condition. This is because the fundamental resonant frequency for 
an open jet is twice as high as the fundamental frequency for a closed wall. The 
onset of resonance may also be detected by the dramatic shift in phase angle be- 
ginning around M=.6 to M=.8. 
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It should be emphasized that the effect of frequency on wind tunnel interference 
is ■ a continuous one except of course at the actual resonant frequencies . For this 
reason, unsteady wind tunnel tests should carefully take into consideration such un- 
steady interference effects even though the reduced frequencies to be tested may lie 
below their resonance values. 


Table 17. vs. M,k and y for an oscillating flap in a wind tunnel 


75% hinge 


n=3. 05556 


Cl- (y=“ - closed wall) 


Cl» ( y=0-open jet) 


0 

.0 

3. 93081 

0.00000 

3.93081 

0.00 

2.45693 

0.00000 

2.45693 

0.00 

0 

.1 

3.63755 

-0.52198 

3.67481 

-8.17 

2.43951 

0.03078 

2.43970 

0.72 

0 

.2 

3.15090 

-0.60877 

3.20917 

-10.94 

2.39288 

0.07723 

2.39413 

1.85 

0 

.3 

2.80157 

-0.45100 

2.83764 

-9.15 

2.32977 

0.14703 

2.33441 

3.61 

2 

.0 

4.01629 

0.00000 

4.01629 

0.00 

2.48840 

0.00000 

2.48840 

0.00 

2 

.1 

3.70295 

-0.56739 

3.74617 

-8.71 

2.47107 

0.02451 

2.47119 

0.57 

2 

. 2 

3.19239 

-0.66681 

3.26128 

-11.80 

2.42506 

0.06541 

2.42594 

1.55 

2 

.3 

2.83538 

-0. 51065 

2.88100 

-10.21 

2.36358 

0.13075 

2.36720 

3.17 

4 

.0 

4. 31034 

0.00000 

4.31034 

0.00 

2.59221 

0.00000 

2.59221 

0.00 

4 

.1 

3.91916 

-0.73544 

3. 98757 

-10. 63 

2.57547 

0.00202 

2.57547 

0.05 

4 

.2 

3. 31937 

-0.87977 

3.43398 

-14.84 

2.53205 

0.02405 

2.53216 

0.54 

4 

. 3 

2.93429 

-0.73451 

3.02483 

-14.05 

2.47706 

0.07086 

2.47808 

1.64 

6 

.0 

4.98528 

0. 00000 

4.98528 

0.00 

2.80743 

0.00000 

2.80743 

0.00 

6 

. 1 

4.35658 

-1.18762 

4.51556 

-15.25 

2.79291 

-0.05360 

2.79343 

-1.10 

6 

.2 

3. 50453 

-1.44000 

3.78884 

-22.34 

2.75829 

-0.08640 

2.75965 

-1.79 

6 

.3 

3.02163 

-1.36395 

3.31521 

-24.30 

2.72364 

-0.08998 

2.72512 

-1.89 

8 

. 0 

6.84313 

0. 00000 

6.84313 

0.00 

3.27336 

0. 00000 

3.27336 

0.00 

8 

.1 

5. 00685 

-2.81467 

5.74377 

-29.34 

3.27198 

-0.22175 

3.27949 

-3.88 

8 

.2 

2.88884 

-3.24022 

4.34102 

-48.28 

3.28232 

-0.45312 

3.31345 

-7.86 

8 

.3 

0.65287 

-2.84171 

2.91574 

-77.06 

3.33763 

-0.74720 

3.42024 

-12.62 

9 

.0 

9.92111 

0.00000 

9.92111 

0.00 

3.79549 

0.00000 

3.79549 

0. 00 

9 

.1 

3.69704 

-5.68263 

6.77941 

-56.95 

3.83945 

-0.51021 

3.87320 

-7.57 

9 

.2 

-0. 75925 

-1.22187 

1.43855- 

121.86 

3.99142 

-1.33569 

4.20898 

-18.50 

9 

.3 

1.80588 

-0.32653 

1.83517 

-10.25 

3.28022 

-4.17219 

5.30726 

-51.83 
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§15. Conclusions 

The niamerical calculation of unsteady airloads in ventilated wind tunnels 
has been extended from our previous work [5] to include airfoils with multiple 
leading and trailing edge flaps. This has been accomplished with further 
developments both in the mathematical theory and in the computer program TWODI. 

The computational foundation has been strengthened by establishing a more 
powerful existence and convergence theory for arbitrary downwashes having finite 
norm. This has led to improved error estimates and to a better understanding of 
the behavior of collocation for solving integral equations with a leading Cauchy 
singularity. Using the method of collocation for discontinuous downwash, our 
theoretical estimates predict errors no better than 0(1/N) where N is the number 
of basis functions. In practice this works out to 1-2% error using collocation 
alone with the maximum number of basis functions (presently twenty-four) . We 
have shown in this work that such errors can be reduced by several orders of 
magnitude using Richardson extrapolation; for essentially the same amount of 
arithmetic, error reduction by factors of 500 have been achieved. Furthermore, 
the method of extrapolation appears to work equally well for steady and unsteady 
flow. To demonstrate this technique, we have presented accurate numerical re- 
sults for an airfoil in a ventilated tunnel with a three quarter chord flap 
oscillating at high frequencies up to and beyond resonance. 

We have also shown that Bland's collocation method can be made much more effi- 
cient by identifying the singularities in the kernel and integrating them separate- 
ly with suitable quadrature rules. In the case of incompressible flow in free air, 
we have demonstrated with TWODI using five basis functions for continuous down- 
washes, that six decimal accuracy can be obtained uniformly for steady and unsteady 
flow with little increase in computer time. 

The applicability of TWODI has been extended by incorporating the kernel (4.13) 
for steady flow in porous wall tunnels. The complete computational solution to this 
problem is presented for unit downwash, and solutions for other downwashes can be 
obtained at will. 

The practical utility of TWODI has been improved by the incorporation of a 
completely new input module. The entire program has been coded in ANSI FORTRAN 
and is available for general use. Complete user instructions and a set of sample 
problems are provided in this report. 
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Several areas of future research appear to be particularly promising at 
the present time. From a practical standpoint, the method of extrapolation 
should be incorporated into the solution process for airfoils with flaps and 
preliminary steps should be undertaken toward automatic accuracy control in 
the future. Additional progress remains to be made through more efficient 
integration of singularities in the kernel, and in the actual computation of 
the kernel as well. The entire problem of boundary conditions at ventilated 
wind tunnel walls remains open to research, and the possibility of using 
nonlinear integral equation methods for transonic flow is an intriguing one. 
Also, the study of solution methods other than collocation, both theoretically 
and computationally, can be expected to be useful in the future. 

While a general and fully satisfactory theory of unsteady wind tunnel 
interference effects does not yet exist, the present work should increase 
our practical computational capability and we hope it will add to the re- 
liability and precision of aerodynamic testing. 

Addendum Added in Proof 

Since the writing of this report, the TWODI program has been extended to 
provide automatic extrapolation and to utilize the improvements arising from 
the reformulation of the Possio kernel described in section 12. These changes 
are reflected only in the user instructions and sample input/output found in 
section 9 and the Appendix. 
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APPENDIX 


This appendix provides example TIMESHARE input and output for three problems. 
Responses by the user are underlined for ease of identification. 

The first problem is entirely defaulted and produces results which agree to 
6 decimals with the exact Sbhngen and Kiissner-Schwarz solutions [5] (Problem one 
is listed to display the initial default values.) The second problem is chosen to 
verify Bland's results [8] for an airfoil oscillating about the 42.5% chord in a 
closed wall tunnel at M=.85. The second problem is not listed and another problem 
is not immediately entered. Computer execution then occurs and the printed output 
follows. 

To demonstrate the modification of old problems into new ones, a new problem 
one is next defined as the old problem two and then edited. The new problem one uti- 
lizes five methods of solution. The first method merely demonstrates checking the 
airfoil polynomials. Methods 2-4, if extrapolated with (11.51), will show that 
of a thin symmetrical airfoil for a flap hinged at the 75% chord in a ventilated 
wind tunnel with n=550/180, y=2.0137, and M=.5 is given by 

CLg = 3.204, (A-1) 

which matches experimental results [68, Table 4] , and is predicted by method 5 
using automatic extrapolation. 
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THIS IS TWODI-III 

FOR IiFFAULT TO INITIAL OR MOST RECENTLY DEFINEli VALUES 
TYPE D FOLLOWED BY CARRIAGE RETURN IF IN TIMESHARE AND ENTER 
AN OTHERWISE BLANK CARD WITH A D IN COLUMN 1 IF IN BATCH. 

IF IN TIMESHARE TYPE HALT TO STOP. 


ARE YOU IN TIMESHARE OR BATCH 


? D 
ENTER 

? ^ 
ENTER 

? JD 
ENTER 

? IL 
ENTER 

? IL 
ENTER 
? D_ 
ENTER 
? JD 
ENTER 
? IL 
ENTER 
? JO. 

ENTER 
? JD 
ENTER 


NUMBER OF LINES PER PACE 
TITLE 

OUTPUT COMBINATION OF FOURIER* SECTION AND WORK 

LIST OF PRESSURE POINTS 

LIST OF MACH NUMBERS 

LIST OF FREOUENCIES 

LIST OF HEIGHT TO CHORD RATIOS 

LIST OF MASS EFFECT VENTILATION COEFFICIENTS 

LIST OF VISCOUS EFFECT VENTILATION COEFFICIENTS 

LIST OF NODES 


? JD 

ENTER NUMBER OF MODE SHAPES 
'i’ JD. 

ENTER MODE SHAPE I 
? JD 

ENTER MODE SHAPE 2 
? JD. 

ENTER MODE SHAPE 3 

? il 

ENTER LIST OF NODE NUMBERS OF HINGES 
? JD 

ENTER NUMBER OF METHODS OF SOLUTION 
T JD^ 

ENTER SOLUTION PARAMETERS FOR METHOD I. 

? iL 

DO YOU WANT THE INPUT DATA LJSTI D 




in I II ini I 
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SAMPLE PROBLEM 

TIMESHARE == T 
LINES PER PAGE - 66 

FOURIER = T 

SECTION = T 

WORK - F 


800000 

-.600000 

10 PRESSURE POINTS 
-.400000 

-.200000 

0 

200000 

.400000 

.600000 

,800000 

1.000000 


1 MACH NUMBERS 
0 

2 REDUCEB FREQUENCIES 
0 

1 .000000 

1 HEIGHT TO CHORIi RATIOS 
INFINITY 

1 WIND TUNNEL MASS EFFECT VENTILATION COEFFICIENTS 

INFINITY 

1 WIND TUNNEL VISCOUS EFFECT VENTILATION COEFFICIENTS 

0 

3 CHORD WISE NODES 

- J . .000000 0 1.000000 


3 MODE SHAPES 


MODE 1 

.S64190 

.564190 

.564190 

MODE 2 

--.564190 

.564190 

1.692569 

MODE 3 

.564190 

-.564190 

2.820948 


0 HINGE LOCATIONS 

1 METHODS OF SOLUTION 
SOLUTION PARAMETERS FOR METHOD 1 

II 3 12 =: 5 13 = 0 R1 ---•= 0 R2 == 0 

:iO YOU WANT TO MAKE CHANGES 
NO 

10 YOU WANT TO ENTER ANOTHER PROBLEM 
YES 
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TO MODIFY OR RETAIN AN OLD PROBLEM? ENTER ITS NUMBER 
OTHERWISE TYPE D FOLLOWED BY CARRIAGE RETURN 
? JD 

ENTER TITLE 

? VERIFICATION OF RESUL TS IN BLAND SJAM J APPL MATH 18(4) 830-848 1970 
ENTER OUTPUT COMBINATION OF FOURIER? SECTION AND WORK 
? SECTION 

ENTER LIST OF PRESSURE POINTS 
? £ 

ENTER LIST OF MACH NUMBERS 
? 1 «85 

ENTER LIST OF FREQUENCIES 
? 3 0 ♦! *2 

ENTER LIST OF HEIGHT TO CHORD RATIOS 
? 1 7.5 

ENTER LIST OF MASS EFFECT VENTILATION COEFFICIENTS 

? El 

ENTER LIST OF VISCOUS EFFECT VENTILATION COEFFICIENTS 
? D. 

ENTER LIST OF NODES 
? 2 -1 1 

ENTER NUMBER OF MODE SHAPES 

i_ 

ENTER MODE SHAPE 1 
? -.85 1>15 

ENTER LIST OF NODE NUMBERS OF HINGES 

? IL 

ENTER NUMBER OF METHODS OF SOLUTION 
? D 

ENTER SOLUTION PARAMETERS FOR METHOD 1 
? D 

DO YOU WANT THE INPUT DATA LISTED 
? NO 

DO YOU WANT TO MAKE CHANGES 

? 

DO YOU WANT TO ENTER ANOTHER PROBLEM 
? NO 
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SAMPLE PROBLEM 


MACH 


0 

SOLUTION BY 


REDUCED 

COLLOCATION. 


P'REQIJENCY 
NUMBER OF BASIS 


0 

FUNCTIONS 


FREE AIR 


FOURIER COEFFICIENTS OF PRESSURE DUE TO DISPLACEMENT MODE 1 


0 0 

0 0 

0 0 

0 0 

0 0 

AERODYNAMIC NORM OF PRESSURE FUNCTION = 0 


FOURIER COEFFICIENTS OF PRESSURE DUE TO DISPLACEMENT MODE 2 


2.00000000 0 

0 0 

0 0 

0 0 

0 0 

AERODYNAMIC NORM OF PRESSURE FUNCTION == 2.00000000 


FOURIER COEFFICIENTS OF PRESSURE DUE TO DISPLACEMENT MODE 3 


6.00000000 0 

-4.00000000 0 

.00000000 0 

.00000000 0 

-.00000000 0 

AERODYNAMIC NORM OF PRESSURE FUNCTION =■• 7.21110255 




SECTION 

AIRLOAD COEFFICIENTS 

MODE 

COEFF 

REAL 

IMAGINARY 

MAGNITUDE 

1 

LIFT 

0 

0 

0 

1 

PITCH 

0 

0 

0 

1 

X CP 

.25000000 

0 

.25000000 

'? 

LIFT 

7.08981540 

0 

7.08981540 


PITCH 

0 

0 

0 


X CP 

.25000000 

-0 

.25000000 

3 

LIFT 

21 .26944621 

0 

21 .26944621 

3 

PITCH 

-3.54490770 

0 

3.54490770 

3 

X CP 

.58333333 

-0 

.58333333 


PHASE ANGL 


-180.00 

0 


LiJCOOOOOO 



w ro 
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MODE 

1 

1 

1 

:l. 

1 

1 

:i 

1 

1 

1 

p 

'1 

2 
9 
9 
2 
2 
2 
3 
3 

3 
3 
3 
3 
3 
3 
3 
3 


PRESSURES 


X 

x/c 

REAL 

IMAGINARY 

MAGNITUDE 

PHASE ANGLE 

-.80000 

♦10000 

0 

0 

0 

0 

-.60000 

.20000 

0 

0 

0 

0 

-.40000 

.30000 

0 

0 

0 

0 

-.20000 

.40000 

0 

0 

0 

0 

0 

.50000 

0 

0 

0 

0 

.20000 

.60000 

0 

0 

0 

0 

.40000 

.70000 

0 

0 

0 

0 

.60000 

.80000 

0 

0 

0 

0 

.80000 

.90000 

0 

0 

0 

0 

1.00000 

1.00000 

0 

0 

0 

0 

-.80000 

. 10000 

13.54055 

0 

13.54055 

0 

-.60000 

.20000 

9.02703 

0 

9.02703 

0 

-.40000 

.30000 

6.89451 

0 

6.89451 

0 

-.20000 

.40000 

5.52791 

0 

5.52791 

0 

0 

.50000 

4.51352 

0 

4.51352 

0 

♦ 20000 

.60000 

3.68527 

0 

3.68527 

0 

.40000 

.70000 

2.95479 

0 

2.95479 

0 

.60000 

.80000 

2.25676 

0 

2.25676 

0 

.80000 

.90000 

1.50451 

0 

1.50451 

0 

1.00000 

1.00000 

0 

0 

0 

0 

-.80000 

.10000 

24.37299 

0 

24.37299 

0 

-.60000 

.20000 

23.47029 

0 

23.47029 

0 

-.40000 

.30000 

23.44134 

0 

23.44134 

0 

-.20000 

.40000 

23.21721 

0 

23.21721 

0 

0 

.50000 

22.56758 

0 

22.56758 

0 

.20000 

.60000 

21.37457 

0 

21 .37457 

0 

.40000 

.70000 

19.50162 

0 

19.50162 

0 

.60000 

.80000 

16.70001 

0 

16.70001 

0 

. 80000 

.90000 

12.33695 

0 

12.33695 

0 

1.00000 

1.00000 

0 

0 

0 

0 
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SAMPLE PROBLEM 


MACH - 0 REDUCED FREQUENCY .1.00000 FREE A.TR 

SOLUT.TON BY COLLOCATION. NUMBER OF BASIS FUNCTIONS ^ 5 


FOURIER COEFFICIENTS OF PRESSURE 

--.39972777 
-.50000591 
. 00000577 
-.000003:1.6 
.00000:1.12 

aerodynamic norm of pressure; 

FI7URIER C0EFFI(.;.[E;NTS of pressure 

. 77942828 
-•.75001865 
-.24998218 
-.00001090 
.00000407 

AERODYNAMIC NORM OF PRESSURE 

FOURIER COEFFICIENTS OF PRESSURE 

3.43720243 
3. 75000:1.26 
- . 4 1 666 7 1 7’ 

- . :l. 6666885 
, 00000178 

AERODYNAMIC NORM OF PRESSURE 


DUE TO DISPLACEMENT MODE 1 
.53943896 
•-.00000177 
.00000160 
-.00000116 
.00000056 

FUNC T I ON = . 83 7 1 2 758 

DUE TO DISPLACEMENT MODE 2 
1 .87833520 
2.00001662 
-.00001614 
.00000834 
-.00000262 

FUNCTION = 2.95985423 

DUE TO DISPLACEMENT MODE 3 
1 . 47722619 
4.00006962 
1 . 99993338 
.00003832 
-.00001229 

FUNCTION = 6.94699559 



SECTION 

AIRLOAD COEF 

FICIENTS 


COEEI- 

REAL 

tMAGINARY 

MAGNITUDE 

PHASE ANGLE 

LIFT 

-1.41699805 

1 .91226 1 32 

2.38004765 

~:l 26.54 

P I rcH 

.44311870 

.00000157 

.4431:1870 

-♦00 

X CP 

.47168992 

. 29917627 

. 55856765 

-32.39 

LIFT 

2.76300130 

6.65852491 

7.20903115 

-67.46 

PITCH 

.66468672 

-1 . 77246858 

1.89300114 

69.44 

X CP 

.63350860 

. 35878941 

.72805424 

-29.53 

LIFT 

12.18456537 

5.23663048 

13.26219937 

-23.26 

PITCH 

-3 . 32335208 

■••3. 54496940 

4.85916424 

133.15 

> CP 

.92154087 

.29326671 

. 96707959 

-17.65 




PRESSURES 


MODE 

X 

X/C 

REAL 

IMAGINARY 

MAGNITUDE 

PHASE ANGLE 

:l 

-*80000 

*10000 

-*67518 

3*65215 

3*71404 

-100*47 

1 

-♦60000 

*20000 

-1*35286 

2*43476 

2*78537 

-119*06 

;t 

-*40000 

*30000 

-1 *72272 

1.85957 

2*53491 

-132*81 

1 

-*20000 

*40000 

-1*93404 

1.49098 

2*44203 

-142*37 

1 

0 

*50000 

-2*03048 

1*21738 

2*36746 

-149.06 

;l 

*20000 

*60000 

-2*02641 

*99399 

2.25706 

-153*87 

:i. 

*40000 

*70000 

-1*92022 

*79696 

2.07904 

-157.46 

1 

*60000 

*80000 

-1.69227 

.60869 

1*79841 

-160*22 

1 

*80000 

*90000 

-1 *27863 

*40579 

1*34148 

-162*39 

1 

1*00000 

1*00000 

0 

0 

0 

0 

2 

“*80000 

*10000 

8*39126 

4*59251 

9*56580 

-28.69 

2 

-*60000 

*20000 

5*05245 

6*67258 

8*36962 

-52*87 

2 

-*40000 

*30000 

3.16940 

7*85409 

8*46947 

-68*02 

2 

- , 20000 

*40000 

1.76728 

8*50845 

8.69005 

-78.27 


0 

*50000 

*63055 

8*75251 

8*77520 

-85*88 

p 

*20000 

*60000 

-.29591 

8*62050 

8*62558 

-91*97 

2 

*40000 

*70000 

-1*00550 

8.09369 

8*15591 

-97*08 

2 

*60000 

*80000 

-1*44499 

7.08436 

7*23022 

-101*53 

2 

*00000 

*90000 

-1 *47486 

5*32471 

5.52519 

-105.48 

2 

1 *00000 

1 * 00000 

0 

0 

0 

0 


-*80000 

* 10000 

7*40125 

-6*78906 
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VERIFICATION OF RESULTS IN BLAND SIAM J APPL MATH 18<4) 830-848 1970 


MACH = .85000 REDUCED FREQUENCY = 0 ETA = 7.5000 

CLOSED TUNNEL WALL 

SOLUTION BY COLLOCATION. NUMBER OF BASIS FUNCTIONS = 5 

SECTION AIRLOAD COEFFICIENTS 

MODE COEFF REAL IMAGINARY MAGNITUDE PHASE ANGLE 

1 LIFT 12.23511824 0 12.23511824 0 

1 PITCH -.03818645 0 .03818645 -180.00 

1 X CP .25624211 -0 .25624211 0 


VERIFICATION OF RESULTS IN BLAND SIAM J APPL MATH 18(4) 830-848 1970 


MACH - .85000 REDUCED FREQUENCY = .10000 ETA •■= 7.5000 

CLOSED TUNNEL WALL 

SOLUTION BY COLLOCATION. NUMBER OF BASIS FUNCTIONS = 5 

SECTION AIRLOAD COEFFICIENTS 

MODE COEFF REAL IMAGINARY MAGNITUDE PHASE ANGLE 

1 LIFT 5.89895953 -5.39528821 7.99417653 42.45 

1 PITCH -.39062630 -.35268436 .52628430 137.92 

1 X CP .26256375 .13106613 .29345878 -26.53 


VERIFICATION OF RESULTS IN BLAND SIAM J APPL MATH 18(4) 830-848 1970 


MACH = .85000 REDUCED FREQUENCY = .20000 ETA = 7.5000 

CLOSED TUNNEL WALL 

SOLUTION BY COLLOCATION. NUMBER OF BASIS FUNCTIONS = 5 

SECTION AIRLOAD COEFFICIENTS 

MODE COEFF REAL IMAGINARY MAGNITUDE PHASE ANGLE 

1 LIFT 5.43100764 -.22007538 5.43546476 2.32 

1 PITCH -.01944215 -.78007066 .78031290 91.43 

1 X CP .24552645 .28708426 .37775734 -49.46 
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DO YOU WANT TO ENTER ANOTHER PROBLEM 

? ^ 

YOU MAY MODIFY OR RETAIN OLD PROBLEM N AS NEW PROBLEM M 
PROVIDED N IS GREATER THAN OR EQUAL TO M 
TO MODIFY OR RETAIN AN OLD PROBLEM^ ENTER ITS NUMBER 
OTHERWISE TYPE D FOLLOWED BY CARRIAGE RETURN 
? 2 

DO YOU WANT THE INPUT DATA LISTED 
? NO 

DO YOU WANT TO MAKE CHANGES 
? YES 

DO YOU WANT TO LINE EDIT 


? YES 

NOW OPEN FOR LINE EDITING. WHEN DONE TYPE END 
? TITLE 
ENTER TITLE 

? COMPARISON WITH TIJDEMAN & SCHIPPERS NLR TR 73018L 1974 
? MACH 

ENTER LIST OF MACH NUMBERS 
? 1 .5 
? FREQUENCY 

ENTER LIST OF FREQUENCIES 

? L_£ 

? HEIGHT 

ENTER LIST OF HEIGHT TO CHORD RATIOS 
? 1 3 .0555556 
T MASS 

ENTER LIST OF MASS EFFECT VENTILATION COEFFICIENTS 
? 1 2.0137 
? NODES 

ENTER LIST OF NODES 
? 3 -1 .5 1 
? MODE 

ENTER NUMBER OF MODE SHAPES 
? 1 

ENTER MODE SHAPE 1 


? 0 0 .5 
? H1N6E" 


ENTER LIST OF NODE NUMBERS OF HINGES 
? 1 2 
? METHOD 


ENTER NUMBER OF METHODS OF SOLUTION 


? 5 

ENTER SOLUTION PARAMETERS FOR METHOD 


? 1 9 0 0 0 


ENTER 

SOLUTION 

PARAMETERS 

FOR 

METHOD 

? 3 3 0 

0 0 




ENTER 

SOLUTION 

PARAMETERS 

FOR 

METHOD 

? 3 6 0 

0 0 




ENTER 

SOLUTION 

PARAMETERS 

FOR 

METHOD 

? 3 12 

0 0 0 




ENTER 

SOLUTION 

PARAMETERS 

FOR 

METHOD 

? 3 3 2 

0 0 





? END 

DO YOU WANT THE INPUT DATA LISTED 
? YES 


1 

2 

3 

4 
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COMPARISON WITH TUriEMAN % SCHIPPERS NLR TR 7301SL 1974 

TIMESHARE ~ T 
LINES PER PAGE - 66 

FOURIER = F 

SECTION - T 

WORK = F 

0 PRESSURE POINTS 

1 MACH NUMBERS 
.500000 

1 REBUCED FREQUENCIES 
0 

1 HEIGHT TO CHORD RATIOS 
3.055556 

1 WIND TUNNEL MASS EFFECT VENTILATION COEFFICIENTS 

2.013700 

1 WIND TUNNEL VISCOUS EFFECT VENTILATION COEFFICIENTS 

0 


- 1 .000000 
MODE 1 

0 


3 CHORDWISE NODES 
.500000 1.000000 

1 MODE SHAPES 

0 .500000 

1 HINGE LOCATIONS 





5 METHODS 

OF SOLUTION 


SOLUTION 

PARAMETERS 

FOR 

METHOD 1 



1 1 -■= 1 

12 = 9 

13 = 

0 R1 = 

0 

R2 = 

SOLUTION 

PARAMETERS 

FOR 

METHOD 2 



1 1 3 

12 = 3 

13 = 

0 R1 = 

0 

K = 

SOLUTION 

PARAMETERS 

FOR 

METHOD 3 



1 1 = 3 

12 = 6 

13 = 

0 R1 = 

0 

R2 = 

SOLUTION 

PARAMETERS 

FOR 

METHOD 4 



11=3 

12 = 12 

13 = 

0 R1 = 

0 

R2 = 

SOLUTION 

PARAMETERS 

FOR 

METHOD 5 



11 = 3 
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TO MAKE CHANGES 
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YOU WANT 
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PRODU H 




? NO 
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***#COMPUTATIONAL CHECK OF MAJOR SUBROUTINES**** 


9*1 CHECKING AFP. 
.2482434168 
.2482434168 


CORRECT VALUES f 
.5941373295 
.5941373295 


THEN CURRENTLY 

-.6544599169 

-.6544599169 


COMPUTED VALUES. 
.5508566527 
.5508566527 


COMPARISON WITH TIJDEMAN S SCHIPPERS NLR TR 73018L 1974 

MACH = .50000 REDUCED FREQUENCY = 0 ETA = 3.0556 

MASS EFFECT (SLOTTED WALL) VENTILATION COEFFICIENT = 2.0137 

VISCOUS EFFECT (POROUS WALL) VENTILATION COEFFICIENT = 0 

SOLUTION BY COLLOCATION. NUMBER OF BASIS FUNCTIONS = 3 


SECTION AIRLOAD COEFFICIENTS 


MODE 

COEFF 

REAL 

IMAGINARY 

MAGNITUDE 

PHASE ANGLE 

1 

LIFT 

2.86621003 

0 

2.86621003 

0 

1 

PITCH 

-.77670622 

0 

.77670622 

-180.00 

1 

X CP 

.79197439 

-0 

.79197439 

0 



COMPARISON 

WITH TIJDEMAN 

8 SCHIPPERS 

NLR TR 730 18L 1974 

MACH = 

.50000 

REDUCED FREQUENCY = 

0 ETA = 

3,0556 

MASS EFFECT 

(SLOTTED WALL) 

VENTILATION 

COEFFICIENT = 

2.0137 

VISCOUS EFFECT 

(POROUS WALL) 

VENTILATION 

COEFFICIENT - 

0 


SOLUTION 

BY COLLOCATION 

!, NUMBER OF 

■ BASIS FUNCTIONS 

= 6 



SECTION AIRLOAD COEFFICIENTS 


MODE 

COEFF 

REAL 

IMAGINARY 

MAGNITUDE PHASE ANGLE 

1 

LIFT 

3.01589975 

0 

3.01589975 

0 

1 

PITCH 

-.74651846 

0 

.74651846 

-180.00 

1 

X CP 

.74505522 

-0 

.74505522 

0 


COMPARISON 

WITH TIJDEMAN 

8 SCHIPPERS 

NLR TR 7301 8L 1974 

MACH = 

.50000 

REDUCED FREQUENCY = 

0 ETA 

3.0556 

MASS EFFECT 

(SLOTTED WALL) 

VENTILATION 

COEFFICIENT - 

2,0137 

VISCOUS EFFECT 

(POROUS WALL) 

VENTILATION 

COEFFICIENT = 

0 


SOLUTION 

BY COLLOCATION 

. NUMBER OF 

BASIS FUNCTIONS 

= 12 


SECTION AIRLOAD COEFFICIENTS 


MODE 

COEFF 

REAL 

IMAGINARY 

MAGNITUDE 

PHASE ANGLE 

1 

LIFT 

3.10514470 

0 

3.10514470 

0 

1 

PITCH 

-.73830865 

0 

.73830865 

-180.00 

1 

X CP 

.72553897 

-0 

.72553897 

0 


COMPARISON 

WITH TIJDEMAN 

8 SCHIPPERS 

NLR TR 730 18L 

1974 

MACH = 

50000 

REDUCED FREQUENCY = 

0 ETA 

= 3.055 

MASS 

EFFECT 

(SLOTTED WALL) 

VENTILATION 

COEFFICIENT = 

2.0137 

VISCOUS EFFECT 

(POROUS WALL) 

VENTILATION 

COEFFICIENT = 

0 


SOLUTION BY EXTRAPOLATION. ORDER = 2 PERIOD 

== 3 



SECTION 

AIRLOAD COEFFICIENTS 


MODE 

COEFF 

REAL 

IMAGINARY 

MAGNITUDE 1 

='HASE ANGLE 

1 

LIFT 

3.20398970 

0 

3.20398970 

0 

1 

PITCH 

-.73468821 

0 

.73468821 

“180.00 

1 

X CP 

.70860835 

-0 

.70860835 

0 

0 YOU WANT 

TO ENTER ANOTHER PROBLEM 




? NO 
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